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I.  INTRODUCTION 

With  the  development  of  the  finite  element  (FE)  method,  and  the 
advancement  of  computer  technology,  it  is  now  possible  to  design  and 
analyze  complex  engineering  systems.   Through  FE  techniques,  a  detailed 
mathematical  model  of  a  complex  structural  dynamic  system  may  be 
developed,  and  the  static  and  dynamic  responses  determined.  The 
'traditional'  procedure  for  conducting  a  finite  element  analysis  (FEA)  of  a 
structure  is  to  first  assemble  the  system  matrices.  These  system  matrices  are 
the  mass,  stiffness,  and  damping  matrices,  which  constitute  the  FE  model  of 
the  complete  structure.  The  next  step  would  be  to  determine  the  system 
responses  using  various  solution  techniques.   This  is  the  most 
computationally  demanding  phase  of  the  FEA  process.   The  results  are  then 
processed  and  interpreted.   As  a  result  of  this  analysis,  the  engineer  may  wish 
to  change  some  aspect  of  the  design  and  perform  the  analysis  again.   Even 
more  useful  would  be  the  use  of  some  type  of  optimization  scheme  to  find 
the  optimum  design  change  while  concurrently  performing  the  FEA. 
However,  if  the  'traditional'  method  of  FEA  is  used,  every  time  a  design 
variable  is  changed,  the  affected  system  matrices  must  be  reassembled  and  the 
entire  solution  phase  must  be  reaccomplished.   Depending  on  the  complexity 
of  the  system,  and  the  number  of  different  design  parameters  which  may  be 
changed,  this  route  would  be  computationally  impractical.   As  a  result,  the 


designer  may  only  be  able  to  iterate  through  the  design  process  a  few  times, 
and  be  left  with  a  design  which  is  less  than  optimum. 

Therefore,  more  efficient  techniques  for  calculating  a  system's 
responses  after  modifications  have  been  made,  must  be  introduced.   One  such 
technique  is  the  use  of  synthesis  procedures  to  obtain  the  modified  system 
responses.   These  synthesis  procedures  involve  the  use  of  the  structure's 
baseline  (pre-modification)  responses,  the  modifications  made  to  the 
structure's  mass,  stiffness,  and  damping  matrices,  and  the  equations  which 
link  the  two  to  obtain  the  modified /synthesized  responses.  The  advantage  of 
this  is  that  the  system  matrix  assembly  and  solution  phase  of  the  FEA  process 
must  only  be  accomplished  once.  A  smaller  set  of  change  matrices  are 
assembled,  and  along  with  the  presynthesized  responses  and  computationally 
efficient  synthesis  equations,  the  synthesized  responses  are  obtained. 

The  two  types  of  sythesis  techniques  to  be  featured  in  this  thesis  are 
frequency  domain  synthesis,  and  time  domain  synthesis.   Frequency  domain 
synthesis  makes  use  of  the  baseline  frequency  response  functions  (FRFs)  to 
calculate  the  new  FRFs  after  the  structure  has  been  modified.   Time  domain 
synthesis  uses  the  impulse  response  functions  (ERFs)  of  the  baseline  structure, 
the  coupling  forces  from  the  modification,  and  the  convolution  integral.  The 
combination  of  these  produces  a  nonstandard  nonhomogeneous  Volterra 
integrodifferential  equation  (VIDE)  of  the  second  kind  which  is  solved  in 
order  to  calculate  transient  responses  of  the  modified  structure.   The  use  of 


these  techniques  greatly  reduce  the  computer  computation  time  and  allow  for 
the  application  of  optimization  techniques  to  complex  engineering  systems. 


II.  FINITE  ELEMENT  FORMULATIONS 

As  mentioned  previously,  the  initial  step  in  the  FEA  process  is  to  create 
a  finite  element  model  of  the  structure  to  be  analyzed.   In  this  study,  different 
finite  element  types  are  used  to  ensure  the  validity  and  universal  applicability 
of  the  synthesis  techniques,  and  to  assist  in  the  synthesis  formulations.  For 
simplicity,  linear  lumped  mass-spring  systems  were  used.   Its  formulation  is 
based  on  the  use  of  Newton's  laws  to  derive  the  equations  of  motion  for  each 
particle  of  mass  [Ref.  1].  These  systems  allowed  for  the  initial  general 
formulation  of  the  synthesis  techniques  and  the  building  of  a  checking  system 
to  ensure  its  accuracy.  In  order  to  illustrate  more  realistic  systems,  systems 
formed  from  beam  elements  and  plate  elements  were  used.   The  beam 
element  formulation  is  based  on  Euler-Bernoulli  beam  bending,  and  the  use 
of  Hermitian  shape  functions  [Refs.  2  &  3].   The  plate  element  formulation  is 
based  on  shear  deformation  which  includes  both  the  transverse  shear  energy 
and  the  bending  energy  [Ref.  3].  One  very  important  anomally  about  the  plate 
shear  deformation  formulation  is  that,  the  unconstrained  stiffness  matrix  (K) 
is  rank  deficient  by  the  number  of  degrees  of  freedom  per  node,  plus  one. 
When  solving  the  eigenvalue  problem  to  obtain  the  system's  natural 
frequencies  (con)  and  mode  shapes,  one  additional  "spurious"  rigid  body  mode 
is  obtained  which  disrupts  all  modal  calculations.   Therefore,  to  avoid  this 
problem,  all  analyses  done  using  the  plate  were  done  on  a  constrained  to 


ground  system  to  eliminate  all  rigid  body  modes.  The  constraints  were  then 
subsequently  synthesized  out  during  the  analysis  in  order  to  obtain  the  truly 
desired  responses. 


III.  FREQUENCY  DOMAIN  STRUCTURAL  SYNTHESIS 

The  following  background  information  on  frequency  domain 
structural  synthesis  is  taken  from  references  [4  &  5].  This  portion  of  the  thesis 
will  outline  the  methodologies  and  formulations  of  the  frequency  synthesis 
technique  and  the  classical  techniques  used  to  check  its  accuracy. 

Frequency  domain  techniques  were  demonstrated  as  early  as  1939  by  G. 
Kron  in  his  book  on  the  tensor  analysis  of  electrical  networks.  Since  then, 
numerous  formulations  and  application  techniques  have  been  used  and 
documented.  The  advantages  in  computational  time,  and  the  flexibility  of 
this  technique,  have  also  been  well  documented  [Ref.  6].   The  solutions 
obtained  through  this  method  are  exact  with  no  approximations. 

As  mentioned  previously,  frequency  domain  structural  synthesis  is  a 
methodology  whereby  changes  can  be  made  to  a  given  structure,  and  its  new 
frequency  response  function  (FRF)  can  be  formulated  through  the  use  of  a 
single  synthesis  equation.  This  is  quite  different  from  the  classical  method  of 
evaluating  a  structural  change  through  the  reconstruction  of  the  basic 
elements  which  define  the  structure,  and  then  using  an  impedance  inversion 
or  modal  superposition  technique  on  the  new  structure.  The  frequency 
synthesis  technique  may  be  divided  into  two  major  classes,  coupling  and 
modification.  Coupling  is  the  joining  of  a  totally  independent  uncoupled 
structure  to  the  given  structure.  Modification  is  the  addition  of  redundant 


load  paths  in  the  given  structure.  Each  of  the  two  classes  of  synthesis  may  be 
either  direct  or  indirect.  For  the  indirect  case,  there  is  the  existence  of  an 
interconnection  impedance  element  between  substructures  in  a  coupling 
operation,  and  degrees  of  freedom  in  a  modification  operation.  In  the  direct 
case,  the  coupling  and  modification  is  done  without  the  impedance  element 
existing,  therefore  making  the  modification  synonymous  to  applying  a 
constraint. 

This  thesis  will  focus  on  the  indirect  modifications  to  a  given  structure. 
It  will  allow  for  the  inclusion  of  a  spring  stiffener  or  a  viscous  damper 
between  two  points,  a  lumped  mass  on  the  structure,  and  implementing  a 
base  excitation  to  the  structure.  A  generalized  computer  program  will  be 
presented,  which  will  allow  for  any  of  the  above  mentioned  changes  to  be 
accomplished  on  some  baseline  structure.  The  program  then  returns  the  FRFs 
for  all  dofs  where  changes  have  occurred  and  any  other  dofs  that  may  be  of 
interest  to  the  user. 

A.    GENERALIZED  FREQUENCY  RESPONSE 

The  initial  step  in  obtaining  the  new  synthesized  FRFs  for  the  modified 
structure,  is  to  have  the  FRFs  for  the  baseline  structure.  The  following 
equations  are  for  a  system  which  is  excited  by  a  force  on  the  structure  or  by  a 
base  excitation.  These  formulations  are  also  used  to  verify  the  accuracy  of  the 
synthesis  after  a  modification  has  been  made. 
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Consider  the  following  two  degree  of  freedom  system: 
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Figure  3.1  Two  Degree  of  Freedom  Mass-Spring-Damper  System 


where: 


mp  =  mass  of  a  plate 

mc  =  mass  of  computer  equipment 

kp  =  isolator  stiffness  between  plate  and  ground 

kc  =  isolator  stiffness  between  plate  and  computer  equipment 

xp(t)  =  plate  displacement  as  a  function  of  time 

xc(t)  =  computer  displacement  as  a  function  of  time 

F(t)  =  excitation  force 


Applying  Newton's  2nd  law,  the  equations  of  motion  for  this  system  may  be 


written  as: 


mp      0 
0      m. 


+ 


c  J  K.     c . 


cn  +  c      — c. 

p  c  c 


■  + 


c   J  K    c  . 


k+kc     -kc 

p  c  c 

-k.        k. 


PI 

101 


(3.1) 


Generalizing  this  formulation  to  an  ndof  (number  of  degrees  of  freedom) 
system  where  the  external  force  may  occur  at  any  dof,  and  using  a  more 
compact  matrix  and  vector  notation,  the  2nd  order  system  of  linear  equations 
for  an  ndof  structure  can  be  written  as: 

[M]{x}  +  [C]{x}  +  [K]{x}  =  {F}  (3.2) 

Assuming  a  harmonic  forcing  function  and  consequently  a  harmonic  steady- 
state  response: 

{F(t)}  =  {F}eiQt  {X(0}  =  {X}^  (3.3) 

where: 

{f}  &  {x}   are  the  amplitudes  of  the  harmonic  forcing  function  and 

response. 
Taking  the  first  and  second  derivatives  of  the  response  (x(t)},  and  substituting 
into  equation  (3.1)  and  simplifying: 

[K  +  JQ.C  -  Q2M]{X}  =  {F}  (3.4) 

[Z(Q)]{X}  =  {F}  (3.5) 

[Z(Q)]  is  known  as  the  impedance  matrix.  Multiplying  both  sides  of  the 
equation  by  [Z(Q)]"1  and  denoting  this  as  frequency  response  function  (FRF) 
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matrix  [H(Q)],  equation  (3.5)  becomes: 

{X}  =  [H(Q)]{F} 

where: 

Hn           H12 

"lxndof 

[H(Q)]  = 

rl21                iT-22 

"2xndof 

■"■ndofxl       "ndof  x  2 

ndof  x  ndof 

(3.6) 


(3.7) 


and  in  elemental  form 


1     E 


(3.8) 


The  FRF  matrix  [H(Q)]  is  both  complex  valued,  and  frequency  dependent.  The 
advantage  of  this  formulation  in  this  work  is  that  it  can  be  used  to  check  the 
frequency  synthesis  method  when  damper  modifications  are  made  to  the 
structure.  These  damper  modifications  represent  non-proportional  damping, 
which  is  not  easily  handled  in  modal  coordinates. 

The  modal  coordinate  formulation  of  the  FRF  matrix  begins  from  the 
same  general  equations  of  motion  presented  in  equation  (3.2),  and  the 
assumptions  of  harmonic  excitation  and  response  in  equation  (3.3).  It  also 
includes  the  assumption  of  a  harmonic  modal  response: 

{*«}  =  fry*  (3.9) 

and  the  linear  modal  transformation: 


{x(t)}  =  [<%«} 


(3.10) 
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where  [O]  is  the  assembles  matrix  of  mass  normalized  mode  shapes.  By 
applying  equations  (3.3),  (3.9),  &  (3.10),  to  equation  (3.2)  and  simplifying: 

[-Q2[M][0]  +  MC][G>]  +  [KW]{q]  =  {F}  (3.11) 

Premultiplying  equation  (3.11)  by  [0]T: 


-Q' 


+  jO. 


2{>, 


CO. 


{<?}  =  {<*}  (3.12) 


where  it  is  recognized  that: 
[0>]T[M][O]  =  [I] 
[0]T[C][0]  =  [2Cco] 
[0]T[K][0)]  =  [co2] 

[<d]t  {f}  =  m 


(3.13) 
(3.14) 
(3.15) 
(3.16) 


£  is  the  modal  damping  factor,  co  is  the  natural  frequency,  and  the  subscript  r 
represents  each  mode.  Rewriting  equation  (3.12)  and  solving  for  the  modal 
displacement  yields: 


w= 


1 


coi2-Q2  +  j2£icoiQ 


{*} 


(3.17) 


Using  equations  (3.10)  &  (3.16),  equation  (3.17)  is  transformed  back  to  physical 
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coordinates: 


M=[o] 


coT  2-Q2+f2Cr(QrQ 


[0]T{p} 


(3.18) 


From  equation  (3.6): 


{H(Q)}  =  [0] 


CQt  z-Q.2  +  j2CIcotQ. 


mT 


(3.19) 


and  any  element  of  the  FRF  matrix  may  be  represented  as: 


nmoQfl 


lihi 


M 


coT2-er+j2Crcorci 


(3.20) 


Although  this  modal  coordinate  formulation  does  not  handle  non- 
proportional  damping,  it  does  provide  for  the  ability  to  sum  over  a  subset  of 
the  complete  modes,  rather  than  summing  over  all  modes  for  the  FRF.  The 
number  of  modes  necessary  to  correctly  capture  the  response  is  dependent  on 
the  frequency  range  of  interest.  The  lower  the  frequency  range  of  interest,  the 
less  number  of  modes  necessary.  The  obvious  advantage  of  this  procedure  is 
that  for  large  dof  systems,  this  summation  may  be  truncated,  and 
computational  time  saved.  However,  the  question  of  how  many  modes  are 
enough  must  be  considered  and  truncation  criteria  must  be  established.  The 
modal  formulation  will  be  the  one  used  for  the  synthesis  techniques. 
Recognizing  the  modal  truncation  issue,  this  thesis  uses  all  modes  in  its  FRF 
calculations. 
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Now  consider  the  following  two  degree  of  freedom  system: 


hVL  r\H 


n\ 


y(t)        Xp(t)  xc(t) 

Figure  3.2  Base  Excited  Two  Degree  of  Freedom  Mass-Spring-Damper  System 


where: 

kb  =  isolator  stiffness  between  plate  and  base 

y(t)  =  base  excitation  displacement 

all  other  variables  are  the  same  as  in  the  system  illustrated  in  Figure 
3.1. 

Applying  Newton's  2nd  law,  the  equations  of  motion  for  this  system  may  be 

written  as: 


mp      0 
0      m. 


+ 


cb  +  cc    -cc 


— c         c 


kb  +  kc    -kc 

-kc      K 


*pl_ 


cb    0 
0     0 


W+ 


kb    o" 


{y} 


0     0 
(3.21) 
Generalizing  this  formulation  to  an  ndof  system  where  the  base  excitation 
may  occur  at  any  dof,  and  using  a  more  compact  matrix  and  vector  notation, 
the  2nd  order  system  of  linear  equations  for  an  ndof  structure  can  be  written 
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as: 

[M]{x}  +  [C]{x}  +  [K]{x}  =  [Cb]{y}  +  [Kb]{y}  (3.22) 

Following  the  same  procedure  as  with  the  FRF  formulation  for  a  force 
excitation  on  the  structure,  and  assuming  a  harmonic  base  excitation 

({y(t)}  =  JYje*01),  the  following  result  is  obtained: 

{X}  =  [K  +  jQC-  n2Mp[Kb  +  £2Cb]{Y}  (3.23) 

Since  the  base  excitation  {yJ  is  the  same  at  all  dofs,  Y  can  be  factored  out  and 
the  result  is: 

{X}  =  [K  +  jQ.C  -  Q2M]_1[Kb  +  *QCb]{l}Y  (3.24) 

and 

{H}  =  [K  +  jQC  -  Q2M]_1[Kb  +  *QCb]{l}  (3.25) 

In  this  instance  the  collection  of  FRFs  is  not  a  matrix  but  rather  a  vector 
where  in  elemental  notation: 

H=^  (3.26) 

Y 

In  some  instances,  the  base  excitation  and  desired  response  may  not  both  be 
displacements.  One  different  combination  could  be  an  interest  in  obtaining 
the  output  displacement  response  due  to  a  given  base  acceleration.  From  the 
assumption  of  a  harmonic  base  excitation  and  response: 

{y(t)}  =  { Y}e^  {x(t)}  =  {X}e**  (3.27a&b) 

{y(t)}  =  {y}^  =  -Q2{Y}e**  {x(t)}  =  |X}e**  =  -Q2{x}^  (3.28a&b) 
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••■  PMa*} 


w~&W 


(3.29a&b) 


The  following  table  uses  the  above  relationships  between  acceleration  and 
displacement  to  illustrate  possible  combinations  and  how  to  handle  each. 


Y 

Y 

X 

M 

M(->/2) 

t 

[H](-G2) 

[H] 

Table  3.1  Response  and  Base  Excitation  Relationships 

The  modal  formulation  for  base  excitation  is  not  needed  since  the 
synthesis  method  for  a  base  excitation  or  force  excitation,  uses  the  previous 
modal  FRF  formulation  exclusively. 

B.         FREQUENCY  DOMAIN  SYNTHESIS  FORMULATION 

Now  consider  the  following  general  ndof  (number  of  degrees  of 
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freedom)  structure,  to  which  structural  modifications  are  to  be  made: 


Figure  3.3  General  NDOF  Structural  System 


The  letters  i,  c,  and  b  represent  the  physical  coordinate  system  for  the 
structure  where: 

i  =  iset  -  the  set  of  internal  dofs  where  no  changes  occur,  but  there  is  an 
interest  in  knowing  FRF  information  about  these  dof. 

c  =  cset  -  the  set  of  dofs  where  changes  have  occurred. 

b  ■  bset  -  the  set  of  dofs  where  the  structure  is  indirectly  connected  to  a 
base  excitation. 
Using  equation  (3.6)  and  the  previously  described  coordinate  system  for 
rearrangement  and  partitioning,  the  structural  system  is  described  in  the 
frequency  domain  at  each  frequency  as: 


V 

<Xc 

xb. 

►  = 

[[HB]j  [HJ 

[Hlb]l 

m 
[4  J 

[Hd] :  [h«] 

[H»] 

[KJ]  ["bcl 

[Hbb]J 

(3.30) 
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Before  modification,  the  force  vector  refers  to  externally  applied  forces 
to  the  structure.  However,  following  modification,  there  exists  coupling 
forces  between  the  cset  and  bset  dofs.  By  definition,  the  iset  dofs  do  not 
experience  these  coupling  forces.  Therefore,  the  partitioned  forces  may  be 
rewritten  as: 


fc=fr+fcz=r-Azcx; 
4=r+fbz=r-Azb(x;-y) 


(3.31a) 
(3.31b) 
(3.31c) 


where: 


f2  -  the  coupling  forces  due  to  structure  modification. 

[AZ]  -  the  impedance  matrix  which  is  equal  to  [AK  +  jQAC  +  Q2AM]. 

x*  -  generalized  (synthesized)  responses  after  modification. 
Substituting  equations  (3.31)  into  equation  (3.30),  and  recognizing  that  the 
presynthesized  responses  {xi  xc  xb}T  are  those  due  to  the  external  forces  only, 
the  following  result  is  obtained: 


*         -s 

* 

*i 

V 

Xc 

► 

=  ■> 

*c 

►  — 

xb. 

xb. 

[Hic]     [Hib] 

[H«]    [H*] 

LP**]    [Hbb]J 


[AZC]        0 
0       [AZJJK 


[Hic]     [H*] 

[H«]    [HA] 

L[Hbc]    [Hbb]J 


[AZC]        0 


0       [AZjJly 


ro' 


(3.32a) 


Denoting  two  additional  coordinate  sets  'e'  and  'z'  where: 


e  =  eset  =  iu  cu  b 


z  =  zset  =  cu  b 
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and  [AE]  is  the  total  impedance  matrix.  Equation  (3.32a)  can  be  rewritten  as: 

W  -W-  [HjAH]{x2}-  +  [H.b][AZk]{y}  (3.32b) 

From  equation  (3.32b)  it  can  be  seen  that  there  are  two  unknown  synthesized 
responses  in  this  one  equation.  Another  independent  equation  can  be 
generated  from  equation  (3.32a)  without  introducing  any  additional 
unknowns.  By  observing  the  bottom  two  rows  of  equation  (3.32a),  the 
following  relationship  is  derived: 

W  =  K}-[H=lAZ]{xJ'+[Ha,lAZl,]{y}  (3.33) 

Solving  for  {x2}*  and  substituting  that  expression  into  equation  (3.32b)  yields 
the  general  expression  for  the  responses  of  the  synthesized  structure. 

M  =  W  -  [HJ[AH]([I  +  H^ASf  {xj  +  [I  +  H^nH^lAZjjy})  +  [Heb][AZb  ]{y} 

(3.34) 
From  the  general  equation  for  frequency  response,  it  is  recognized  that: 

W-W)  K}  =  [H«]{fr}  (3.35a&b) 

Equation  (34)  is  therefore  rewritten  as: 

{xj*  =  [H.]{r }  -  [H  J[AH]([I  +  HJ&F[Hj[fr}  +  [I  +  H^AHnH^lAZjM)  +  [Heb][AZb]{y} 

(3.36) 
As  can  be  seen  from  equation  (3.36),  the  synthesized  response  is  now  in  terms 
of  the  known  FRFs  of  the  presynthesized  structure,  known  impedance  due  to 
modifications,  and  known  force  and  base  excitations. 
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Now  that  the  complete  general  frequency  response  synthesis  equation 
has  been  derived,  consider  the  case  where  there  is  no  base  excitation  applied 
to  the  structure.  By  letting  {y}  =  0,  the  general  synthesis  equation  becomes: 

W  =  ([Hee]  -  [H  J[AH][I  +  H^AiFfH  J){r  }  (337) 

Since  no  base  excitation  exists,  the  impedance  between  the  structure  and  base 
(AZb),  is  reduced  to  a  V  dof  to  ground  change.  Therefore  the  bset  dofs  are 
really  additions  to  the cset:  z=cub=cuc=c and  AZ^  =  AZC.  Also 

recognizing  that  {xe}*  =[Hee]*{feext},  the  synthesized  FRF  matrix  at  each 
frequency  can  be  written  as: 

[Hj*  =  [H  J  -  [HjfAZjl  +  HccAZc;T[Hce]  (3.38) 

Now  consider  the  case  where  there  are  no  external  forces  applied  to  the 
structure.  By  letting  fff)  =  0,  the  general  synthesis  equation  becomes: 

{"«}'  =  ([HejAZb]-[Hez][AHlI  +  Hz7AH]-'[Hzb][AZb]){y}  (3.39) 

Factoring  out  the  base  excitation  magnitude  y  and  applying  equation  (3.6)  yet 
again: 

{He}*  =  ([H.b][AZb]-[He2][AH][l  +  HzzAZ]-,[H!blAZb]){l}       .         (3.40) 

In  elemental  notation,  the  FRF  vector  represents  the  following  at  each 
frequency: 

{H.}-={x.}'I  (3.41) 
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Equations  (3.38)  and  (3.40)  show  how  the  new  synthesized  FRFs  are 
obtained  from  the  presynthesized  FRFs  and  the  impedance  caused  by 
structural  modifications.  These  synthesis  equations  are  exact,  and  there  are  no 
limitations  to  the  modification  values.  The  change  may  even  be  negative  to 
represent  the  removal  of  mass  or  removal  of  an  interconnection  element 
from  the  structure. 

C         FREQUENCY  DOMAIN  SYNTHESIS  COMPUTER  CODE 

The  computer  language  used  for  the  frequency  domain  synthesis  and 
all  other  computer  coding  in  this  thesis  is  MATLAB  V.4.2c.l.  The  goal  in  the 
formulation  of  the  frequency  domain  synthesis  computer  programs  is  to 
perform  the  synthesis  for  an  indirect  modification  operation  to  a  given 
structure.  The  program  allows  for  the  inclusion  of  a  spring  stiffener  or  a 
viscous  damper  between  two  dofs,  a  lumped  mass  addition  on  the  structure, 
and  implementing  a  base  excitation  to  the  structure.  The  program  then 
returns  the  FRFs  for  all  dofs  where  changes  have  occurred  and  any  other  dofs 
that  may  be  of  interest  to  the  user.  Appendix  A  shows  the  computer  codes 
used  to  perform  the  frequency  domain  synthesis,  and  perform  a  comparative 
analysis  of  the  synthesis  versus  classical  methods.  The  term  classical  is 
synonymous  to  the  traditional  method  of  reformulating  the  elemental 
matrices  following  a  modification. 
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The  programs  presented  in  Appendix  A  can  be  segregated  into  three 
categories.  The  first  is  the  calculation  of  the  synthesized  FRFs  due  to  an 
external  force  excitation  on  the  structure,  and  their  comparison  to  the 
classically  calculated  FRFs.  The  second  is  the  calculation  of  the  synthesized 
FRFs  due  to  a  base  excitation  on  the  structure,  and  their  comparison  to  the 
classically  calculated  FRFs.  Both  of  these  programs  use  the  impedance 
inversion  method  to  calculate  the  presynthesized  FRFs,  and  the  classical  FRFs 
following  modification.  The  third  category  is  the  calculation  of  the 
synthesized  FRFs  where  the  different  aspects  of  the  synthesis  technique  have 
been  modularized  into  MATLAB  functions.  This  modularization  assists  the 
process  in  running   more  efficiently,  and  allows  for  the  universal  use  of  the 
frequency  synthesis  technique.  This  is  crucial  since  the  ultimate  goal  of  the 
technique  is  to  be  used  efficiently  in  another  process  (i.e.  optimization).  Since 
efficiency  is  of  concern,  the  modal  method  of  calculating  the  presynthesized 
FRFs  is  used. 

The  following  is  a  diagram  of  the  flowpath  of  the  modularized 
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frequency  synthesis  programs: 


LoadFRFs 
changcm 

»  ' 

Load  Structure 
changejn 

-> 

Calculate  FRFs 

modal.m 
frfmodal.m 

-i 

Modifications 
changejn 

Form  Change 
Matrices 

deltajn 

—» 

Perform  Synthesis 
frfsynthjn 

FRF  Selection 
frf  siftjn 


Figure  3.4  Flowpath  of  the  Modularized  Frequency  Synthesis  Programs 

A  more  detailed  description  of  each  function  is  contained  within  Appendix 
A. 


D. 


ILLUSTRATIVE  EXAMPLES 


The  purpose  of  this  section  is  to  illustrate  the  use  of  the  frequency 
domain  synthesis  technique  and  validate  its  accuracy.  The  computational 
efficiency  of  this  technique  has  already  been  well  documented  [Ref.  6], 
therefore  a  time  comparison  is  not  done.  The  three  structural  systems 
previously  described,  mass-spring-damper,  beam  and  plate,  are  used  to 
accomplish  this  goal. 
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1.         Mass-Spring-Damper  System 

Consider  the  following  mass-spring-damper  system: 


~24 


." lilt '  k=  1000  lb/in 

k                 j                         1  k^  =  250  lb/in 

-f/ft—IIIH^IIII-T-IIIM  «—             m  =  0.259  lb-s2/in 

m                     ,    !  F£(t)            c30  =  21b-s/in 


C30 


u 


X(t) 

Figure  3.5  Mass-Spring-Damper  System  Experiencing  Force  Excitation 

The  solid  elements  represent  the  original  structure  and  the  dashed  elements 
represent  modifications.  The  subscript  'i'  in  the  excitation  force  represents  the 
dof  where  the  excitation  is  applied.  Although  not  shown  in  the  figure,  the 
original  structure  is  proportionally  damped  using  [C]  =  a[K].  The 
modifications  are  arbitrarily  selected  and  the  addition  of  the  damper  element 
represents  nonproportional  damping.  The  following  is  the  resultant 
synthesized  and  classically  determined  FRF  H32(£2)  which  represents  the 
response  magnitude  at  dof  3  due  to  a  force  at  dof  2. 
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Frequency  Response  Function  (H32) 


4  6 

Frequency  (Hz) 


Figure  3.6  Force  Excitation  Mass-Spring-Damper  System  FRF 

From  the  graphs  produced,  it  can  been  seen  that  the  plots  are  identical, 
proving  that  the  synthesis  technique  is  exact,  and  that  the  computer  coding 
for  the  arbitrary  changes  is  correct.  These  results  were  produced  using  the 
fsynstr.m  program. 
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2.  Beam  System 

Now  consider  the  following  beam  system: 

^_^_^_^_^_^_^_^_^_^  ^  =  1000  lb/in 

[ [ | ;| •   '.  >;   r     » » | ;'  c130  =  20  lb-s/in 

p  /       i  f  mn  =  .5  lb-s2/in 

4  ».  Kj,  mn     i_     ,  c130  >  . 

Ty(t)  «i  T  ^    _yC0J 

St* 

Figure  3.7  Beam  System  Experiencing  Base  Excitation 

where  the  beam  parameters  are  : 

Length:  L  =  5  ft  Width:  w  =  3  in  Height:  h  =  4  in 

Young's  Modulus:  E  =  10e6  psi 

density:  p  =  2.53e-4  lbf-secA2/inA4 
The  free-free  beam  which  is  discretized  into  10  beam  elements,  11  nodes,  and 
22  dofs,  represents  the  original  structure.  The  modifications  are  represented 
by  the  dashed  lines.  Just  as  in  the  previous  example,  proportional  damping  is 
used  to  form  the  original  [C]  matrix  and  all  changes  to  the  beam  are  arbitrary. 
The  following  is  the  resultant  synthesized  and  classically  determined  FRF 
H13(&)  which  represents  the  response  magnitude  at  dof  13  due  to  the  base 
excitations. 
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Frequency  Response  Function  (H13) 


400  600 

Frequency  (Hz) 


800 


1000 


Figure  3.8  Base  Excitation  Beam  System  FRF 

Just  as  before,  it  can  be  seen  that  the  two  plots  are  identical.  These  results  were 
produced  using  the  fsynbase.m  program. 

3.         Mass-Plate  System 

Up  until  now  the  structures  use  to  demonstrate  the  frequency  synthesis 
technique  were  relatively  simple  with  relatively  small  number  of  dofs.  The 
true  test  of  the  technique  and  the  computer  code  is  in  their  applicability  to 
relatively  large  structures.  Therefore,  consider  the  following  computer-plate 
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system: 


Ak,,  =  le4  lb/in 
Ac^  =  1  lb-s/in 
Akc  =  500  lb/in 
Ac  =1 


&*y 


2 


Figure  3.9  Plate-Mass  System  Experiencing  Base  Excitation 


Thickness:  t  =  1  in 
Poisson's  Ratio:  v  =  0.3 


where  the  computer  and  plate  parameters  are  : 
Plate:       Width:  w  =  5  ft  Depth:  d  =  5ft 

Young's  Modulus:  E  =  30e6  psi 

Density:  p  =  7.35e-4  lbf-secA2/inA4 
Computer:        Weight:  W  =  100  lbs 

In  this  case,  the  original  structure  is  represented  by  the  free-free  plate 
with  the  attached  computer  located  off  center  and  above  the  plate.  The  plate  is 
discretized  into  25  plate  elements,  36  nodes,  and  108  dofs.  The  computer  is 
modeled  as  a  single  lumped  mass  in  the  translational  direction  only,  and 
increases  the  total  number  of  nodes  and  dofs  to  37  and  109  respectively.  The 
modifications  are  represented  by  the  dashed  lines  and  are  comprised  of 
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changes  to  the  isolators  between  the  base  and  the  plate.  Just  as  in  the  previous 
example,  proportional  damping  is  used  to  form  the  original  [C]  matrix,  and 
any  changes  to  the  system  are  arbitrary.  However,  since  the  manipulation  of 
this  system  in  an  optimization  routine  is  the  ultimate  goal  of  this  research, 
the  isolator  elements  which  will  act  as  optimization  variables  are  chosen  for 
modification.  Therefore,  modifications  were  also  performed  on  the  isolators 
between  the  plate  and  the  computer.  The  modularized  program  which  uses 
the  modal  method  of  calculating  the  original  FRF  will  be  used  in  this 
instance.  However,  the  classical  method  will  also  be  used  as  a  check  system. 
The  following  is  the  resultant  synthesized  and  classically  determined  FRF 
H109(£2)  which  represents  the  response  magnitude  at  dof  109  (the  computer) 
due  to  the  base  excitations. 
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Frequency  Response  Function  (H109) 
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Figure  3.10  Free-Free  Plate-Mass  System  FRF  (Modal  Synthesis) 

From  the  plots  it  is  easily  observed  that  the  synthesis  method  and  the 
classical  method  do  not  yield  the  same  results.  After  successfully  performing 
countless  other  frequency  synthesis  analysis  on  all  three  structures,  and 
successfully  re-running  this  same  analysis  where  the  original  FRF  is  obtained 
using  impedance  inversion  methods  vice  modal  methods  (Fig.  3.11  shows), 
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Figure  3.11.  Free-Free  Plate-Mass  System  FRF  (Impedance  Synthesis) 

it  is  concluded  that  the  problem  lies  in  the  use  of  the  modal  method  of 
calculating  the  original  FRF.  After  observing  the  mode  shapes  and  natural 
frequencies  of  the  original  free-free  plate-mass  system,  the  existence  of  an 
additional  spurious  rigid  body  mode  is  discovered.  The  existence  of  this  rigid 
body  mode  vice  a  flexible  mode  makes  the  FRF  calculations  incorrect.  After 
unsuccessful  attempts  at  trying  to  replace  the  rigid  body  modes  naturally 
generated  by  this  finite  element  formulation,  with  rigid  body  modes  created 
using  the  Graham-Schmidt  method  [Ref.  2],  the  decision  to  constrain  the 
plate-mass  and  then  remove  the  constraints  using  synthesis  is  reached.  What 
this  accomplishes  is  to  eliminate  the  rigid  body  modes  entirely  and  allows  for 
the  proper  calculation  of  the  original  FRF.  However,  there  is  a  computational 
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price  to  pay  with  the  addition  of  four  more  load  paths  to  contend  with  in  the 
calculations.  This  is  however,  a  very  minute  price  to  pay  when  compared  to 
using  the  impedance  inversion  method. 

The  following  is  the  resultant  modally-synthesized  and  classically 
determined  FRF  H109(Q),  when  the  original  plate-mass  structure  is  exactly  as 
before,  but  with  1000  lb /in  spring-to-ground  elements  located  at  all  four 
corners.  In  the  process  of  the  analysis,  these  elements  are  synthesized  out  of 
the  structure.  Other  modification  values  are  still  the  same  as  before. 
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Figure  3.12  Corners-to-Ground  Plate-Mass  System  FRF  (Modal  Synthesis) 

It  can  be  seen  from  these  plots  that  eliminating  the  rigid  body  modes  in  the 
original  structure,  and  then  synthesizing  the  constraint  elements  out  is  an 
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accurate  means  of  avoiding  the  problem  with  the  rigid  body  modes.  Another 
verification  of  the  ability  to  synthesize  out  elements  is  that  fact  that  these 
plots  also  match  the  plots  in  Figure  3.11  (free-free  plate-mass  system);  as  well 
they  should.  This  will  therefore  be  the  method  used  in  order  to  handle  the 
dynamic  analysis  of  the  plate-mass  system. 
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IV.  STATIC  DISPLACEMENT  SYNTHESIS 

When  using  optimization  procedures  for  the  design  of  a  system,  there 
are  usually  various  aspects  of  the  design  problem  which  are  at  odds  with  one 
another.  For  example,  in  the  design  of  a  support  beam,  the  aspect  to  be 
optimized  may  be  the  weight  of  the  beam.  One  method  to  reduce  weight 
would  be  to  reduce  the  beam's  cross  sectional  area.  However,  since  this  is  a 
support  beam  and  will  therefore  be  carrying  some  sort  of  load,  a  reduction  in 
cross  section  results  in  an  increase  in  stress.  Since  there  are  material 
limitations  on  the  amount  of  stress  the  beam  can  experience,  the  reduction  in 
cross  sectional  area  is  constrained  by  the  stress  limit. 

For  a  shock  and  vibration  isolation  system,  one  aspect  of  the  design 
which  should  be  considered  is  the  static  sag  or  displacement  of  the  system. 
The  static  sag  could  be  used  in  the  optimization  problem  as  a  constraint 
which  limits  the  amount  the  isolator  stiffness  may  be  modified.  Since  the 
solution  for  the  static  response  of  a  structure  involves  the  inverse  of  the 
stiffness  matrix  ([K]"1),  for  large  structures,  it  would  not  be  computational 
efficient  to  perform  this  operation  every  time  the  optimization  variables  are 
changed.  Therefore,  some  sort  of  synthesis  technique  must  be  employed  to 
alleviate  the  need  to  calculate  [K]"1  directly.  This  portion  of  the  thesis  explores 
the  methodologies  and  formulations  for  the  use  of  frequency  domain 
synthesis  in  the  calculation  of  static  displacement.  From  henceforth,  this 
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technique  will  be  referred  to  as  static  displacement  synthesis.  Classical 
techniques  such  as  Guyan  model  reduction  are  used  to  check  the  accuracy  of 
these  formulations. 

A.        GENERALIZED  STATIC  DISPLACEMENT 

From  Newton's  2nd  law,  the  static  equations  for  an  ndof  system  is: 

{F}  =  [K]{x}  (4.1) 

By  defining  the  force  vector  as  the  weight  of  the  structure,  and  assuming  that 
there  are  no  external  forces  on  the  structure,  the  displacement  vector  {x} 
represents  the  static  displacements  of  the  structure  due  to  its  own  weight. 
Therefore,  the  static  displacements  of  the  structure  due  to  its  own  weight 
equals: 

Mstat  =  [K]-:{F}=  KT'EMHg}  (4.2) 

where  {g}  represents  a  vector  where  the  gravitational  constant  appears  at  all 
dofs  that  are  affected  by  gravity  (i.e.  vertical  translational  dofs). 

Equation  (4.2)  returns  the  static  displacements  at  all  dofs  of  the 
structure.  The  size  of  the  stiffness  matrix  of  which  the  inverse  is  taken  is  ndof 
by  ndof.  Since  the  repetitive  calculation  of  [K]"1  is  unsatisfactory  due  to  the 
possibility  of  its  large  size,  and  a  relatively  small  subset  of  the  static 
displacements  of  the  structure  are  desired,  the  use  of  model  reduction  and 
static  condensation  schemes  are  investigated. 
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B.         GUYAN  REDUCTION  /  STATIC  CONDENSATION 

The  following  concepts  and  equations  for  model  reduction  are  taken 
from  references  [2  &  7].  The  Guyan  reduction  defines  a  transformation  matrix 
[T]  which  can  operate  on  a  subset  of  the  structure's  displacements  and  return 
the  full  ndof  set  of  displacements  for  the  structure. 


Xa 


=  [T]{xa}  (4.3) 


where: 


"KqoK^ 


(4.4) 


a  =  aset  -  those  dofs  arbitrarily  selected  as  the  set  of  active  dofs. 

o  ■  oset  -  the  remainder  of  the  dofs  omitted  from  the  aset. 
Substituting  equation  (4.3)  into  equation  (4.1),  and  premultiplying  both  sides 
of  the  equation  by  [T]T  yields: 

{F}  =  [K]{xa}  (4.5) 

where: 

[K]  =  [T]T[K][T]  (4.6a) 

{F}  =  [T]T{F}  (4.7a) 

Therefore,  the  static  displacements  of  the  desired  dofs  equal: 

{xa}  =  [K]"1{F}  (4.8) 
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It  can  be  seen  from  equation  (4.8)  that  the  main  cost  of  obtaining  the  desired 
static  responses  of  the  system  is  now  the  inverse  of  the  reduced  stiffness 

matrix  [K]  which  has  size  of  aset  by  aset  vice  the  inverse  of  the  full  stiffness 
matrix  which  has  size  ndof  by  ndof.  There  is  a  definite  improvement  in 
computational  efficiency,  but  the  issue  of  forming  the  transformation  matrix 
[T],  which  contains  the  inverse  of  the  sub-matrix  [KJ  still  exists.  For  a  large 
structure,  the  size  of  the  aset  will  usually  be  substantially  less  than  the  size  of 
the  oset,  thereby  making  [K^]"1  a  time  demanding  operation.  In  an 
optimization  procedure,  [K^]'1  can  be  performed  prior  to  the  iterations  begin, 
and  passed  into  the  iteration  portion  of  the  optimization  program.    Even 
though  [K^]"1  must  only  be  performed  once,  the  formation  of  [T]  is  still 
undesirable  due  to  the  possibility  that  FRF  data  is  readily  accessible  and  the  [K] 
matrix  is  unavailable.  Therefore,  a  method  of  obtaining  the  reduced  stiffness 
matrix  and  reduced  force  vector  from  FRF  data  is  desired. 

C         STATIC  DISPLACEMENT  SYNTHESIS  FORMULATION 

During  the  derivation  of  the  frequency  domain  synthesis,  the  eset 
coordinate  set  was  defined  as  iucub.  This  means  that  the  eset  includes  the  set 
of  all  dofs  that  the  user  is  interested  in,  dofs  where  modifications  occur,  and 
dofs  where  a  base  excitation  is  attached.  In  other  words,  when  looking  at  the 
structure  as  a  whole,  eset  is  synonymous  to  the  set  of  active  dofs  or  aset. 
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1.         Formulation  of  Reduced  Stiffness  Matrix  [K*] 

Using  the  aset/oset  coordinate  system  previously  described  but  with 
the  eset  equaling  aset,  the  expanded  version  of  equation  (4.6a)  is  written  as: 

[K]  =  [Kee-KeoK00-1Koe]  (4.6b) 

Now  using  the  same  eset/oset  coordinate  system  and  the  fact  that  [H]  =  [Z]"1, 
the  following  relationship  is  written: 


[H.]  !   [Hj 


.[Hoe] !  [HJJJZJ  !  [Zj_ 


I?JJJ?-J 


[I] 

6 


o 
Pi 


(4.9) 


Carrying  out  the  matrix  multiplication  of  equation  (4.9),  the  first  row  yields 
the  following  two  relationships: 

HeeZee+HeoZoe=I  HeeZeo  +  HTOZoo=0  (4.10a&b) 

Solving  for  H^  in  equation  (4.10b)  and  substituting  it  into  equation  (4.10a) 
yields: 


HeeZee  ~  HJLJL^    Z^  -  I 


(4.11) 


Solving  for  Hce  yields: 

[Hee]  =  [Zee-ZeoZ00-1Zoe]"1  (4.12) 

Recognizing  that  Z  =  K  +jQC  -  Q2M,  and  applying  the  static  condition  that 
Q=0: 


[HeeiO^^-KjAr 


(4.13) 
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From  equations  (4.6b)  &  (4.13),  it  is  determined  that: 

[K]  =  [Hee(0)]_1  (4.14a) 

Therefore,  since  the  reduced  stiffness  matrix  can  be  determined  from  an  FRF 
matrix,  when  a  structure  is  modified,  the  new  reduced  stiffness  matrix  can  be 
determined  from  a  synthesized  FRF  matrix. 

[r]  =  [H^(0)]_1  (4.14b) 


2.         Formulation  of  Reduced  Force  Vector  {F*} 

Recognizing  that  {F}  =  [M]{l}g  the  expanded  version  of  equation  (4.7a) 
in  eset/oset  partitioning  is: 

{F}  =  [MM  -  KAX  J  M^  -  KJK^MjWg  (4-Tb) 

Another  very  important  concept  of  Guyan  reduction  is  the  reduced  mass 
matrix. 


[M]  =  [T]T[M][T] 


I 


ee 


By  comparing  equations  (4.7b)  and  (4.15b),  it  can  be  seen  that: 

{F}  =  [M]{le}g 
if  and  only  if 


ee 
"Koo     Koe 


{U=m{ie}= V 


(4.15a) 
(4.15b) 


(4.16) 


(4-17) 
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where  {le}  is  an  eset  length  vector  of  all  ones,  and  {10}  is  an  oset  length  vector 
of  ones  and  zeros,  with  ones  at  the  translation  dofs  and  zeros  at  the  rotations. 
In  effect,  what  equation  (4.17)  is  stating  is  that  for  every  row  of  the 
transformation  matrix,  the  sum  of  its  columns  is  exactly  one  or  zero.  In 
general,  this  occurrence  is  not  true.  It  is  very  dependent  on  which  dofs  are 
chosen  as  the  eset.  Therefore,  it  must  be  determined  which  dofs  are  allowable 
choices  for  the  eset  in  order  for  this  formulation  to  work. 

Consider  the  general  stiffness  matrix  of  a  restrained  structure: 

KR  =  Ku  +  KG  (4.18a) 

where  Ku  represents  the  stiffness  matrix  of  the  unconstrained  structure  and 
KG  represents  the  grounded  stiffness  matrix.  Equation  (4.18a  )  can  be  rewritten 
in  eset/oset  partitioning  format  as: 


[K*]  = 


_L 


*2 


+ 


KJi_.o. 

L  °  I  IK~]J 


(4.18b) 


The  zeros  in  the  grounded  matrix  are  due  to  the  fact  that  grounds  produce  no 
off  diagonal  information.  Now  a  displacement  vector  is  chosen  which 
satisfies  the  static  equations  of  equilibrium  where  only  forces  appear  in  the 
active  or  eset  dofs,  and  which  produces  no  strain  energy  internal  to  the 
restrained  structure. 


[^]{g=[Ki{;:}+[Kip:}= 


[<\       0 

o    K\ 


10 


(4.19) 
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From  equation  (4.19)  it  can  be  seen  that: 


(4.20a) 


The  vector  which  satisfies  equation  (4.19)  is  already  chosen  by  default.  The 
vector  -J    e  >  equals  the  previously  defined  <  e  >.  Therefore,  {*Fe}  =  {le}  and 

Now  consider  the  following  portion  of  equation  (4.19): 


kl.9_: 


w. 


L  o  ■[K^jKrioj 

Repartitioning  this  equation  into  translational  and  rotational  dofs: 


"[k£]t      o         o 

0        [K^]R        0 

0         o      [i&l 


0 


0 


0 


0 

'{%W 

f{F.G}xl 

0 

{*.}. 

►  =  i 

{?}.  . 

0 

Wl 

0 

Kk 

.W.. 

0 

where: 


Wt' 

Wt 

{*.}. 

w. . 

w, 

{UT 

.{*.}«. 

I  W.J 

(4.19) 


(4.21) 


(4.22) 


Case  I  -  Taking  the  third  row  of  equation  (4.21)  yields: 

[KyTMT=[KyT{l„}T  =  {0}  (4.23) 

Since  [K^,]  is  a  diagonal  matrix,  the  only  way  equation  (4.22)  is  true  is  if 

[K^]t  =[0].  Therefore,  this  proves  that  there  cannot  be  any  translational 
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grounds  in  the  oset.  Consequently,  this  means  that  ail  translational  grounds 
must  be  in  the  eset. 

Casell  -  Taking  the  fourth  row  of  equation  (4.21)  yields: 

Equation  (4.24)  is  always  true  regardless  of  [K^l  .  Therefore,  this  shows  that 

there  are  no  restrictions  on  rotational  grounds  being  in  the  oset.  Thus  far  it 
has  been  determined  what  must  be  in  the  eset,  all  translational  grounds.  It 
must  now  be  determined  what  else  is  allowed  in  the  eset. 

Consider  the  zero  internal  strain  energy  requirement: 


I 
I 

t 
I 


K" 


K" 


I*. 


[01 

101 


(4.20b) 


Taking  the  first  row  of  equation  (4.20b)  following  matrix  multiplication 
yields: 

K]w=-K]« 

Repartitioning  this  equation  into  translational  and  rotational  dofs: 


(4.25a) 


i 

i 

TTj_ 
i 

.ft! 


RR. 


eo 


*  IT. 
RT 


JR. 
RR. 


(4.25b) 


Casein  -  if  eset  includes  every  translational  dof  and  no.  rotational  dofs, 
equation  (4.25b)  becomes: 


0         0 


'WtI_ 
.ft}.. 


o  K] 

0         0 


[{0} 


(4.26) 
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After  matrix  multiplication,  equation  (4.26)  requires  that: 

KLMr  =  W  <4-27> 

which  means  that  for  every  row,  the  sum  of  the  translational  dof  columns  for 
an  unrestrained  structure  equals  zero.  This  condition  is  always  true. 
Therefore,  this  shows  that  all  translational  dofs  are  allowed  to  be  included  in 
the  eset. 

Case  IV  -  if  eset  includes  every  translational  dof  and  one  rotational  dof, 
equation  (4.25b)  becomes: 


RT 


Tl_ 


"o 

K] 

TR 

0 

K\ 

RR_ 

'Mr 


(4.28) 


After  matrix  multiplication  and  the  fact  that  [K^]    {le}T  =  {0},  equation  (4.28) 
requires  that: 

[k«UUr={0}  (4-29) 

Since  there  is  only  one  rotation  in  the  eset,  the  dimensions  of  [K^l     are  the 

#  of  eset  translations  by  1.  Consequently,  in  order  for  equation  (4.29)  to  be 
true,  [K^.]     must  always  equal  zero.  By  the  definition  that  this  case  includes  a 

rotational  dof  value  in  the  eset,  [K^]     *[0].  Therefore,  one  rotational  dof 

may  not  be  included  in  the  eset. 

Case  V  -  if  eset  includes  every  translational  dof  and  greater  than  one 
rotational  dof,  the  same  requirement  (equation  4.29)  as  stated  in  Case  IV 

exists.  However,  for  this  case,  [K^.1     does  not  have  to  identically  equal  zero, 
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but  the  sum  of  the  columns  for  every  row  of  [K^]     must  equal  zero.  In 

general,  this  sum  is  not  zero.  Therefore,  it  is  concluded  that  no  rotational  dofs 
can  be  included  in  the  eset. 

Case  VI  -  if  eset  includes  a  subset  of  all  translational  dofs,  equation 
(4.25b)  becomes: 


Kir    0 


M 


TI- 


KI        [K"l 
0  0 


m 


w 


o     oJIMr 

After  matrix  multiplication  and  rearrangement,  equation  (4.30)  requires  that: 

[KSL{Ut+KUUt  =  {<>}  (4.31) 

Equation  (4.31)  is  the  same  as  if  summing  all  translational  columns  of  an 
unrestrained  structure  for  every  dof  in  the  eset.  From  Case  in  (equation  4.27), 
it  is  seen  that  this  sum  is  zero  and  equation  (4.31)  is  satisfied.  Therefore,  this 
shows  that  any  subset  of  all  the  translational  dofs  may  be  included  in  the  eset. 
As  a  result  of  Cases  I-VT,  the  following  requirements  on  the  choices  of 
eset  dofs  for  static  displacement  synthesis  are  summarized: 

(1)  eset  must  include  all  translational  grounded  dofs 

(2)  structure  may  be  grounded  anywhere 

(3)  any  subset  of,  or  all  translational  dofs  may  be  included  in  the  eset 

(4)  no.  rotational  dofs  are  allowed  in  the  eset 

For  the  work  to  be  presented  in  this  thesis,  the  requirements  necessary  for  the 
reduced  force  vector  to  be  determined  from  the  reduced  mass  matrix  are  met. 
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Therefore,  the  task  is  now  to  derive  a  methodology  for  extracting  [M]  from 
FRF  data. 

3.         Formulation  of  Reduced  Mass  Matrix  [M* ] 

Since  [HJO)]  =  [K]"\  it  follows  that: 

[HJQ)]  =  [k-Q2MTl  (4.32) 

Taking  two  derivatives  of  equation  (4.32)  with  respect  to  Q  yields: 

d  He%(fl)  =  [K  -  Q2MTl  [2QM]  ^S^  +  ( [K  -  Q'M]"1  [2M]  +  ^^- [20M]\k  -  Q2MT' 
dQ  dQ      V  dQ  J 

(4.33) 
Applying  the  static  condition  of  Q=0  to  equation  (4.33),  and  solving  for  [M]: 


I[K]^(0)- 

2  dQ2 


[M]  =  t[KJ        ~vw[Kj  (4.34a) 


or  rewritten  in  terms  of  the  synthesis  process: 

[Ml  =  I[h;(0)]-'  ^fflfH^O)]-1  (4.34b) 

From  equations  (4.34a&b),  the  unknown  second  derivative  of  the  eset  FRF  at 

zero  frequency  is  present.  The  calculation  of  this  value  will  be  handled  using 

the  following  forward  finite  differencing  scheme  with  error  order  0(AQ2): 

d2H^(0)  =  -H;(Q3)-h4H^(Q2)-5H;(Q1)  +  2H;(0) 
dQ2  (AQ)2 

The  accuracy  of  this  calculation  is  very  sensitive  to  the  size  of  AQ.  The  proper 

choice  of  AQ  is  dependent  on  the  relative  magnitudes  of  the  mass  and 
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stiffness  matrices.  It  was  found  for  this  work  that  AQ.  =  0.1  provided  a  very 


d2H^(0) 

da2 


accurate  estimation  of  — 7^ — .  However,  this  value  is  not  universal  and  the 


proper  selection  of  AQ.  should  be  based  on  individual  structures. 

From  the  previously  discussed  requirements  on  the  conversion  of  [M] 

to  {F},  it  can  be  seen  that  static  displacement  synthesis  is  not  arbitrary  in  the 
selection  of  which  dofs  may  be  included  in  the  eset.  However,  this  restriction 
does  not  apply  to  which  dofs  can  be  chosen  to  be  modified,  but  only  to 
whether  static  displacement  information  can  be  obtained  about  these  dofs.  For 
example,  modifications  to  rotational  dofs  may  be  made  to  a  structure  and  the 
synthesized  FRF  obtained.  As  a  result  of  this,  rotational  dofs  exist  in  the  eset, 
which  is  not  allowed.  The  solution  to  this  dilemma  is  to  re-synthesize  the 
structure  with  no  changes  and  only  translational  dofs  in  the  iset.  What  this 
does  is  remove  the  rotational  dofs  from  the  eset  and  condense  the  FRFs  prior 
to  using  them  to  obtain  the  reduced  stiffness  matrix  and  force  vector. 
Therefore,  the  static  displacement  method  is  really  not  restrictive  in  its 
application,  but  rather  in  what  information  it  can  provide. 

D.        STATIC  DISPLACEMENT  SYNTHESIS  COMPUTER  CODE 

The  goal  in  the  formulation  of  the  static  displacement  synthesis 
computer  programs  is  to  first  perform  the  frequency  domain  synthesis  for  an 
indirect  modification  operation  to  a  given  structure.  The  program  allows  for 
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the  same  type  of  modifications  to  the  structure  that  were  presented  in  the 

frequency  domain  synthesis  chapter  of  this  thesis.  The  program  then  uses  the 

synthesized  FRFs  and  returns  the  static  displacements  for  all  dofs  designated 

in  the  eset. 

The  reason  why  the  frequency  domain  synthesis  is  performed  here 

seperately  vice  using  the  FRFs  obtatined  from  doing  a  frequency  domain 

synthesis  problem  is  two  fold.  First  of  all,  there  is  no  damping  effect  for  a 

static  problem.  Therefore,  the  synthesized  FRFs  must  be  generated  with  zero 

damping.  Secondly,  if  we  want  to  determine  the  static  displacement  for  a  base 

excitation  problem,  the  base  excitation  form  of  frequency  domain  synthesis  is 

not  what  is  desired.  The  correct  form  of  the  frequency  domain  synthesis  is  the 

force  excitation  synthesis  equation,  with  the  base  excitation  interconnection 

elements  considered  as  springs  to  ground.  Furthermore,  an  entire  spectrum 

of  FRF  information  is  not  necessary  for  static  displacement  calculations. 

Primarily  only  the  static  condition  of  Q=0  is  needed.  However,  in  this 

implementation,  the  synthesized  FRFs  at  the  first  four  frequencies  are  needed 

d2H*  (0) 
for  the  finite  differencing  calculation  of ^ — . 

Appendix  B  shows  the  computer  codes  used  to  perform  the  static 
displacement  synthesis,  and  perform  a  comparative  analysis  of  the  synthesis 
versus  classical  Guyan  reduction  methods.  The  programs  presented  in 
Appendix  B  can  be  seperated  into  two  categories.  The  first  category  is  a 
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program  which  uses  the  classical  Guyan  reduction  techniques  to  calculate 
static  displacement.  The  second  category  is  the  calculation  of  the  synthesized 
static  displacements  where  the  different  aspects  of  the  synthesis  technique 
have  been  modularized  into  MATLAB  functions.  The  following  is  a  diagram 
of  the  flowpath  of  the  modularized  static  displacement  synthesis  programs: 


LoadFRFs 

change.m 

* 

' 

Load  Structure 

change.m 

-> 

Calculate  FRFs 

mod  aim 
frfmodal.m 

— ; 

Modifications 

> 
statchange.m 

— S 

Form  Impedance 

statdelta.m 

— 9 

Perform  Synthesis 

stats  ynth.m 

Figure  4.1  Flowpath  of  the  Modularized  Static  Displacement  Synthesis 

Programs 

It  can  be  seen  that  the  flowpaths  for  static  displacement  and  frequency  domain 
synthesis  are  very  similar.  This  is  done  intentionally  so  that  all  synthesis 
techniques  can  be  driven  by  the  same  modifications  to  a  structure.  A  more 
detailed  description  of  each  function  is  contained  within  Appendix  B. 

E         ILLUSTRATIVE  EXAMPLES 

The  purpose  of  this  section  is  to  illustrate  the  use  of  the  static 
displacement  synthesis  technique  and  validate  its  accuracy.  The  mass-plate 
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structural  system  used  in  the  frequency  domain  synthesis  example,  along 
with  the  changes  that  were  made  to  the  structure,  is  also  used  here  to 
illustrate  static  displacement  synthesis.  The  following  tables  are  comparisons 
of  the  reduced  matrices,  vectors,  and  static  displacements  calculated  using 
frequency  synthesis  and  Guyan  reduction. 
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[K*]  =  [H^(0)f 


.Oe+04* 

1.1518 

-0.0000 

-0.0000 

0.2276 

-0.3794 

0.0000 

-0.0000 

1.0000 

-0.0000 

-0.0000 

-0.0000 

0.0000 

-0.0000 

-0.0000 

1.0000 

-0.0000 

-0.0000 

0.0000 

0.2276 

-0.0000 

-0.0000 

1.3415 

-0.5691 

0.0000 

-0.3794 

-0.0000 

-0.0000 

-0.5691 

1.4985 

-0.5500 

0.0000 

-0.0000 

0.0000 

0.0000 

-0.5500 

0.5500 

Table  4.1a  Static  Displacement  Synthesis  Reduced  Stiffness  Matrix 


[K]  = 

[T]T[K][T] 

1.  Oe+04  * 

1.1518 

0.0000 

-0.0000 

0.2276 

-0.3794 

0 

0.0000 

1.0000 

-0.0000 

0.0000 

-0.0000 

0 

-0.0000 

-0.0000 

1.0000 

-0.0000 

0.0000 

0 

0.2276 

0.0000 

0.0000 

1.3415 

-0.5691 

0 

-0.3794  -0.0000  -0.0000  -0.5691  1.4985  -0.5500 

0  0  0  0  -0.5500  0.5500 


Table  4.1b  Guyan  Reduction  Reduced  Stiffness  Matrix 
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d2H^(0)  =  -H^(^)  +  4H^(Q,)-5H^(Q1)  +  2H^(0) 
dQ2  (AQ)2 


1.0e-07  * 


0.0671 

0.0294 

0.0294 

0.0271 

0.1094 

0.1470 

0.0294 

0.0588 

0.0147 

0.0294 

0.0593 

0.0593 

0.0294 

0.0147 

0.0588 

0.0294 

0.0593 

0.0593 

0.0271 

0.0294 

0.0294 

0.0774 

0.1414 

0.1978 

0.1094 

0.0593 

0.0593 

0.1414 

0.3567 

0.5049 

0.1470 

0.0593 

0.0593 

0.1978 

0.5049 

0.8242 

Table  4.2a  Static  Displacement  Synthesis  2nd  Derivative  Matrix 


t^(0)=[K][2M][K] 


dQ2 


1.0e-07  * 

0.0671 

0.0294 

0.0294 

0.0271 

0.1094 

0.1470 

0.0294 

0.0588 

0.0147 

0.0294 

0.0593 

0.0593 

0.0294 

0.0147 

0.0588 

0.0294 

0.0593 

0.0593 

0.0271 

0.0294 

0.0294 

0.0774 

0.1414 

0.1978 

0.1094 

0.0593 

0.0593 

0.1414 

0.3568 

0.5049 

0.1470 

0.0593 

0.0593 

0.1978 

0.5049 

0.8242 

Table  4.2b  Guyan  Reduction  2nd  Derivative  Matrix 
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J~2L    ~     J 

dQ2     L    " 

'J 

0.1928 

0.0903 

0.0903 

-0.0492 

0.0422 

-0.0000 

0.0903 

0.2940 

0.0735 

0.0620 

0.1417 

-0.0000 

0.0903 

0.0735 

0.2940 

0.0620 

0.1417 

-0.0000 

-0.0492 

0.0620 

0.0620 

0.1537 

-0.0095 

-0.0000 

0.0422 

0.1417 

0.1417 

-0.0095 

0.4215 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

0.2588 

Table  4.3a  Static  Displacement  Synthesis  Reduced  Mass  Matrix 


[M]  = 

[T]T[M][T] 

0.1928 

0.0903 

0.0903 

-0.0492 

0.0422 

0 

0.0903 

0.2940 

0.0735 

0.0620 

0.1417 

0 

0.0903 

0.0735 

0.2940 

0.0620 

0.1417 

0 

-0.0492 

0.0620 

0.0620 

0.1537 

-0.0095 

0 

0.0422 

0.1417 

0.1417 

-0.0095 

0.4215 

0 

0 

0 

-0 

0 

0 

0.2588 

Table  4.3b  Guyan  Reduction  Reduced  Mass  Matrix 
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Synthesis:     {F*}  =  [M*]{le}g  Guyan:  {F}  =  [T]T{F} 

-141.6000  -141.6012 

-255.6080  -255.6102 

-255.6080  -255.6102 

-84.5960  -84.5966 

-285.0197  -285.0226 

-99.9984  -100.0000 


Table  4.4  Reduced  Force  Vectors 


Synthesis:     fc}^  =£KT{F*}  Guyan:     {xe}stat  =  [K]"1{F} 

-0.0296  -0.0296 

-0.0256  -0.0256 

-0.0256  -0.0256 

-0.0316  -0.0316 

-0.0714  -0.0714 

-0.0895  -0.0895 


Table  4.5  Static  Displacements  (in.) 
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From  Tables  4.1  -  4.5,  it  can  be  seen  that  the  static  displacement  synthesis 
calculations  match  those  of  the  Guyan  reduction  at  every  step  of  the  process 
on  its  way  to  calculate  static  displacement. 
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V.  TIME  DOMAIN  STRUCTURAL  SYNTHESIS 

The  following  background  information  on  time  domain  structural 
synthesis  is  taken  from  reference  [8].  This  portion  of  the  thesis  will  outline 
the  methodologies  and  formulations  of  the  time  synthesis  technique  and  the 
classical  techniques  used  to  check  its  accuracy.  There  will  also  be  a  time 
comparison  study  done  in  order  to  quantify  the  computational  efficiency  of 
the  synthesis  technique. 

As  mentioned  previously,  time  domain  structural  synthesis  is  a 
methodology  whereby  changes  can  be  made  to  a  given  structure,  and  its  new 
transient  response  determined  as  the  solution  of  a  nonstandard 
nonhomogeneous  Volterra  integrodifferential  equation  (VIDE)  of  the  second 
kind.  This  integral  equation  is  the  result  of  the  combination  of  the  impulse 
response  functions  (IRFs)  of  the  baseline  structure,  the  reaction  forces 
generated  from  the  modification  to  the  structure,  and  the  total  dynamic 
response  of  the  original  system  in  terms  of  the  convolution  integral.  This 
greatly  differs  from  the  classical  method  of  evaluating  a  structural  change.  In 
the  classical  method,  the  basic  elements  which  define  the  structure  are 
reconstructed,  and  then  a  modal  superposition,  convolution  integral  or  direct 
integration  technique  is  used  on  the  new  structure  to  solve  for  the  response. 

The  motivation  for  the  development  of  the  time  domain  structural 
synthesis  technique  stems  from  the  advantages  and  flexibility  that  is  enjoyed 
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through  the  use  of  the  frequency  domain  structural  synthesis  technique.  The 
time  domain  synthesis  can  be  divided  into  the  same  two  classes  of 
modification  and  coupling,  direct  or  indirect.  Some  of  the  advantages  and 
flexibility  shared  by  the  two  methods  are:  1)  the  formation  of  an  exact 
solution,  with  the  ability  to  handle  damping  modifications;  2)  arbitrary  model 
reduction  in  the  sense  that  only  information  about  those  dofs  which  are 
needed  is  retained;  3)  the  solution  phase  of  the  finite  element  analysis  is  only 
done  once,  and  the  new  synthesis  model  elemental  or  system  matrices  are 
never  formed. 

Just  as  with  the  frequency  domain  synthesis,  this  thesis  will  only  focus 
on  indirect  modifications  to  a  given  structure.  It  will  also  allow  for  the 
inclusion  of  a  spring  stiffener  or  a  viscous  damper  between  two  points,  a 
lumped  mass  on  the  structure,  and  implementing  a  base  excitation  to  the 
structure.  A  generalized  computer  program  will  be  presented,  which  will 
allow  for  any  of  the  above  mentioned  changes  to  be  accomplished  on  some 
baseline  structure.  The  program  then  returns  the  transient  response  for  all 
dofs  where  changes  have  occurred  and  any  other  dofs  that  may  be  of  interest 
to  the  user. 

A.        GENERALIZED  TRANSIENT  RESPONSE 

The  following  equations  are  for  determining  the  total  dynamic 
response  of  a  system  which  experiences  some  sort  of  excitation.  These 
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formulations  are  taken  from  references  [1  &  2],  and  are  used  to  verify  the 
accuracy  of  the  synthesis  after  a  modification  has  been  made. 


1.         Convolution  Integral  Method 

The  convolution  integral  equation  for  determining  the  total  dynamic 
response  of  a  system  is  as  follows: 

{x(t)}  =  {x(t)}h  +  {[hCt  -  T)]{F(T)}dT  (5.1) 

where: 

{x(t)}h  -  homogeneous  solution  which  contains  the  constants  of 
integration 

h(t)  -  matrix  of  impulse  response  functions  (IRFs) 

F(t)  -  excitation  force 

x  -  dummy  time  variable 

The  IRF  matrix  is  determined  modally  using  the  equation: 


[h(t)]  =  [<X>] 


i    -cr<y . 

— e     r    r  C1 


COd  t 

r 


sina)d  t 

r 


[O] 


(5.2) 


where  [O]  =  the  assembled  matrix  of  mass  normalized  mode  shapes.  Any 
element  of  the  IRF  matrix  may  be  represented  as: 

i    -cr<y 

o        X      I    « 

i    t 

r 


Vt)=  ILm—e  "r   rsin*>drt 


(5.3) 
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This  is  a  relatively  efficient  means  of  calculating  a  structure's  dynamic 
response.  However,  this  method  is  based  on  the  principal  of  superposition 
which  is  valid  for  linear  systems  only.  Also  the  modal  method  of  obtaining 
the  IRFs  does  not  readily  handle  non-proportional  damping.  Therefore,  if  a 
damping  change  to  the  structure  is  made,  or  if  a  non-linear  interconnection 
element  is  used  (time  domain  synthesis  with  non-linear  elements  is  possible, 
but  will  not  be  explored  in  this  study),  the  response  cannot  be  determined 
using  the  convolution  integral  method. 

2.         Direct  Integration  Method 

The  direct  integration  method  of  calculating  a  system's  dynamic 
response  is  able  to  handle  the  shortcomings  experienced  by  the  convolution 
integral  method.  However,  the  price  is  paid  in  the  computational  time 
necessary  for  this  solution.  For  the  work  illustrated  in  this  thesis,  MATLAB's 
ODE45.m  function  is  used  for  this  calculation  [Ref.  9].  In  order  to  use  this 
function,  the  2nd  order  ordinary  differential  equations  (ODEs)  of  motion: 

[M]{x}  +  [C]{x}  +  [K]{x}  =  {F}  (3.2) 

must  be  rewritten  as  a  system  of  1st  order  ODEs.  The  following  is  the  general 
matrix  equation  which  converts  any  linear  system's  equations  of  motion  to  a 


60 


1st  order  system. 


fell 


'[I]     o" 

-i 

[B- 

"0    41]" 

JK]     [C]_ 

I{Z,}1" 

(5.4) 


where: 

{zj  =  {x}  {z,}  =  {x} 

{z2}  =  {*}  {z2}  =  {x} 

ODE45.m  solves  for  the  vector  [{zj  {Zj}]7  which  are  the  system's  response 
displacements  and  velocities  using  an  adaptive  time  stepping  procedure. 

B.         TIME  DOMAIN  SYNTHESIS  FORMULATION 

Consider  again  the  exact  same  general  ndof  structure  which  was 
introduced  in  Chapter  3  /  Section  B,  and  to  which  structural  modifications  are 
to  be  made: 


Figure  5.1  General  NDOF  Structural  System 
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Using  the  convolution  integral,  equation  (5.1),  and  the  previously 
defined  sets,  iset,  cset,  bset,  zset,  eset  for  partitioning,  the  total  dynamic 
response  of  the  system  can  be  written  as: 


xi 

xi 

Xc 

►  =  < 

Xc 

►  + 

xb. 

xb. 

h 

I' 

Jo 


[h.a-^lthjt-^J^a-T)] 

"[hd(t-T)I;lhcc(t-T)];[hcb(t-T)] 
[Xi(t-r)] ;  [hJ(V-V)J;r[hb;(t"T)] 


F^T) 
FC(T) 
Fb(T) 


dr 


(5.5) 


Before  modification,  the  force  vector  refers  to  externally  applied  forces 
to  the  structure.  However,  following  modification,  there  exists  coupling 
forces  between  the  cset  and  bset  dofs.  By  definition,  the  iset  dofs  do  not 
experience  these  coupling  forces.  Therefore,  the  partitioned  forces  may  be 
rewritten  as: 


Fi=F;ext 

f] = f-* + f; = Ff  -  (amcX; + accX; + akc<) 

Fb  =  Fr  +  Fb*  =  FT  -  [ACb(x;  -  y)  +  AKb(X;  -  y)] 


(5.6c) 
(5.6b) 
(5.6c) 


where: 


f  -  the  coupling  forces  due  to  structure  modification. 

[AM],  [AC],  [AK]  -  the  modification  matrices 

x*,  x* ,  x* ,  -  generalized  (synthesized)  responses  after  modification. 
Substituting  equations  (5.6)  into  equation  (5.5),  and  recognizing  that  the 
presynthesized  responses  {Xj  xc  xb}T  are  those  due  to  the  external  forces  only, 
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the  following  result  is  obtained: 


bj 


ki  ft[K<t-*>];  [hjt-x)]" 

- k  -     K(t-TJ];  [h4(t-r)j 


VL 


[AMC]    0 
0       0 


■fx;l   [[acc]     o  If  x;  1  [[akc]     o|x;| 
.l°J  L  °    [ACffc-yJ  L  °     [^b]Jlx;-yj 


dT 


(5.7) 


Equation  (5.7)  is  the  nonstandard  nonhomogeneous  VIDE  of  the  second  kind. 
In  order  to  obtain  the  solution  for  this  equation,  the  integral  is  discretized 
using  a  trapezoidal  quadrature  rule.  Also  by  assuming  that  there  are  no 

external  excitations  to  the  structure,  {Feext}  =  0,  equation  (5.7)  becomes: 


LxbJ 


[Aic] 


[AJ 


[Aib] 


[A*] 
[AJ 


F     ^ 


v,  "riCrK,, 


(5.8a) 


or 


The  elemental  form  of  the  quadrature  matrix  is 


(5.8b) 


N= 


t(N), 


2 
At 


0 


0 


At 


0 
0 


M,  f(h,)„      o     o 


0 


(5.9) 


.fW.  MnL  MnL  -  f K)o. 

where  the  subscript  n  represents  the  number  of  time  steps,  and  the  coupling 
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forces  equal: 


0 


AC 


AC„ 


AKC 


AKv 


[AMC]     0" 
0        0 


0 


[ACC]        0 
0        [ACb]^ 


.xb-yj 


[AKC]        0 
0        [AKb]_ 


(5.10a) 


(5.10b) 


(5.10c) 


Equation  (5.8b)  can  be  rearranged  and  placed  in  the  general  form  of 
[A]{x*}  =  {b}  where  {x*}  can  be  solved  for  directly  using  various  methods. 
However,  due  to  the  possible  large  size  of  [A],  it  is  advantageous  to  leave 
equation  (5.8b)  as  is,  and  solve  for  {x*}  by  iteration.  The  iteration  procedure  is 
as  follows: 

1)  assume  the  force  vector  {F2} 

2)  calculate  the  response  vector  {xe}* 

3)  use  {xj*  to  calculate  {xe}\  {xe}\and  a  new  {Fz} 

4)  use  the  new  {Fz}  and  calculate  a  new  {xj* 

5)  repeat  steps  3  &  4  until  {xj*new  -  {xe}*old  <  convergence  tolerance 

Although  the  time  domain  structural  synthesis  is  exact  in  its 
formulation,  approximations  have  now  been  introduced  through  the 
solution  techniques  of  the  trapezoidal  quadrature,  and  the  iteration  method. 
Not  only  is  the  convergence  important,  but  the  stability  of  this  solution  is  also 
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of  concern.  The  stability  and  convergence  is  dependent  on  the  norm  of  Aez 
and  the  magnitude  of  the  modifications  AM,  AC,  AK.  The  norm  is  defined  as 
the  largest  singular  value  of  the  matrix  and  in  this  case  is  driven  by  the  value 
of  At.  The  smaller  the  value  of  At,  the  better  the  stability  and  faster 
convergence  is  reached.  The  larger  the  value  of  the  modifications,  the  greater 
the  possibility  for  instability  in  the  solution.  Therefore,  due  to  the  solution 
technique,  there  are  now  limitations  on  the  modification  values.  Inherent  to 
and  dependent  on  the  computer  system  used,  there  will  be  limitations  on  the 
amount  of  time  which  can  be  covered  by  this  analysis  due  to  memory 
capabilities. 

C         TIME  DOMAIN  SYNTHESIS  COMPUTER  CODE 

The  goal  in  the  formulation  of  the  time  domain  synthesis  computer 
programs  is  to  perform  the  synthesis  for  an  indirect  modification  operation  to 
a  given  structure,  allowing  the  same  type  of  modifications  presented 
previously  in  the  frequency  domain  synthesis  chapter  of  the  thesis.  The 
program  then  returns  the  transient  response  for  all  dofs  where  changes  have 
occurred  and  any  other  dofs  that  may  be  of  interest  to  the  user.  Appendix  C 
shows  the  computer  codes  used  to  perform  the  time  domain  synthesis,  and 
perform  a  comparative  analysis  of  the  synthesis  versus  classical  methods. 

The  programs  presented  in  Appendix  C  can  be  segregated  into  two 
categories.  The  first  is  the  calculation  of  the  synthesized  transient  responses 
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due  to  a  base  excitation  on  the  structure,  and  their  comparison  to  the 
classically  calculated  responses.  One  of  the  programs  in  this  category  uses  the 
convolution  integral  method  in  order  to  check  the  synthesis,  while  the  other 
uses  direct  integration.  The  second  category  is  the  calculation  of  the 
synthesized  transient  responses  where  the  different  aspects  of  the  synthesis 
technique  have  been  modularized  into  MATLAB  functions.  The  following  is 
a  diagram  of  the  flowpath  of  the  modularized  time  domain  synthesis 
programs: 


LoadRFs 
change.m 


Load  Structure 
changcm 


Calculate  IRFs 

modal  .m 
impmodaLm 


Modifications 
changcm 


Form  Change 
Matrices 

delta.m 


Input  Base 
Excitation 

fBlastForcingjn 


Build  Quadrature 
Matrix 

buildA.m 


Perform  Synthesis 


timesynth.m 


Response 
Selection 

time  siftjn 


Figure  5.2  Flowpath  of  the  Modularized  Time  Synthesis  Programs 


A  more  detailed  description  of  each  function  is  contained  within  Appendix  C. 
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D.         ILLUSTRATIVE  EXAMPLES 

The  purpose  of  this  section  is  to  illustrate  the  use  of  the  time  domain 
synthesis  technique  and  validate  its  accuracy.  A  computational  time 
comparison  will  be  done  between  the  synthesis  and  classical  methods  in  order 
to  assess  the  computational  efficiency  of  this  technique.  The  three  structural 
systems  previously  described,  mass-spring-damper,  beam  and  mass-plate,  are 
used  to  accomplish  this  goal. 

1.         Mass-Spring-Damper  System 

Consider  the  following  mass-spring-damper  system: 


m2                 k^ 

k  =  1000  lb/in 

k\:             : 

m  =  0.259  lb-s2/in 

Ak^  =  250  lb/in 

—WHSHM-rllll-* 

Akb  =  250  lb/in 

i  — '  —  <       m 

Am2  =  .02  lb-s2/in 

cb!-'-!     •}-*:  kb 

Ac.  =  2  lb-s/in 

y(t)^    -y-      '; 

y(t)  =  Y0(e-100t  -  e300t; 
Y0=l 

Figure  5.3  Mass-Spring-Damper  System  Experiencing  Base  Excitation 

The  solid  elements  represent  the  original  structure  and  the  dashed  elements 
represent  modifications.  Although  not  shown  in  the  figure,  the  original 
structure  is  proportionally  damped  using  [C]  =  ct[K].  The  modifications  shown 
are  arbitrarily  selected. 
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a.         Spring  Modification  Only 

The  following  figure  is  the  resultant  synthesized  and  classically 
determined  transient  response  x2(t)  which  represents  the  response  at  dof  2 
due  to  the  base  excitation.  For  this  example  neither  the  mass  nor  the  damper 
modification  was  included. 


0.3 


0.2  - 


0.1 


c 
e 

§   o 

8 

JS 

Q_ 


-0.1  - 


-0.2  - 


-0.3 


Transient  Time  Response  at  dof  2 


■  ■          ■-  ■  i 

■  ■■ 

r     —       ->-^ 

i                j 

L. 

Classical 

Synthesis 

0.05 


0.1 
time  (sec) 


0.15 


0.2 


Figure  5.4  Transient  Response  Mass-Spring-Damper  System 

From  the  graph  produced,  it  can  been  seen  that  the  responses  are  identical, 
proving  that  the  synthesis  technique  is  exact,  and  that  the  computer  coding 
for  the  arbitrary  changes  is  correct.  These  results  were  produced  using  the 
tsynconv.m  program. 
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h  Spring  &  Mass  Modifications 

The  previous  example  is  now  repeated  but  with  the  mass  change 
at  dof  2,  Am2  =  .02  lb-s2/in  included.  The  following  figure  shows  the  results  of 
this  additional  modification. 
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Figure  5.5  Transient  Response  Mass-Spring-Damper  System  w/ Spring-Mass 

Modification 

From  the  graph  produced,  it  can  been  seen  that  the  mass  modification  had  no 
effect  on  the  accuracy  of  the  solution.  However,  the  addition  of  the  mass 
modification  did  have  a  large  impact  on  the  computational  time  required  for 
the  synthesis  method.  This  will  be  discussed  in  more  detail  later. 
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c  Spring  &  Damper  Modifications 

The  spring  modification  only  example  is  now  repeated  but  with 
a  damper  change  at  the  base  Acb  =  2  lb-s/in  included.  The  following  figure 
shows  the  results  of  this  additional  modification. 
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Figure  5.6  Transient  Response  Mass-Spring-Damper  System  w/ Spring- 
Damper  Modification 

From  the  graph  produced,  it  can  been  seen  that  the  responses  are  identical 
with  respect  to  the  resolution  of  the  plot  graphics.  Also,  to  further  validate 
the  correctness  of  the  computer  algorithms,  the  effect  of  the  added  damper 
can  be  seen  when  comparing  Fig.  5.4  &  Fig.  5.6.  Since  the  non-proportional 
damper  change  was  made  to  the  structure,  the  direct  integration  method  was 
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used  to  obtain  the  classical  solution.  Although  there  was  no  real  effect  in  the 
accuracy  of  the  solution,  the  effect  in  computational  time  is  profound  for  the 
classical  method.  These  results  were  produced  using  the  tsynode45.m 
program. 

d.         Spring ,  Mass,  &  Damper  Modifications 
The  example  is  now  repeated  with  all  modifications  included. 
The  following  figure  shows  the  results  of  these  modification. 
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Figure  5.7  Transient  Response  Mass-Spring-Damper  System  w/ Spring-Mass- 
Damper  Modification 
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Just  as  in  all  the  other  cases,  it  can  been  seen  that  the  responses  are  extremely 
close  to  identical.  These  results  were  also  produced  using  the  tsynode45.m 
program. 

2.         Beam  System 

Now  consider  the  following  beam  system: 
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Figure  5.8  Beam  System  Experiencing  Base  Excitation 

The  beam  parameters  and  finite  element  modeling  of  the  above  structure  is 
identical  to  that  used  in  the  frequency  domain  synthesis  chapter  of  the  thesis. 
The  modifications  shown  are  arbitrarily  selected. 

a.         Base  Spring  Modification  Only 

The  following  figure  is  the  resultant  synthesized  and  classically 
determined  transient  response  xn(t)  which  represents  the  response  at  dof  11 
due  to  the  base  excitations.  As  can  be  seen  from  Fig.  5.7,  this  is  not  a  dof  where 
change  has  occurred  or  where  the  base  excitation  is  attached,  but  just  a  dof 
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where  there  is  interest  in  the  response.  For  this  example  the  damper 
modification  Ac130  was  not  included. 
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Figure  5.9  Transient  Response  Beam  System 

From  the  graph  produced,  it  can  been  seen  that  the  responses  are  not 
identical,  but  follow  each  other  with  relatively  little  error.  The  difference 
between  the  two  methods  can  be  attributed  to  the  numerical  solution  of  the 
synthesis  method.  As  the  time  step  for  analysis  is  decreased,  the  numerical 
solution  approaches  the  exact  solution.  However,  the  use  of  smaller  time 
steps  is  restricted  by  computer  memory  limitations. 
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h.         Damper  Modification 

The  previous  example  is  now  repeated  but  with  the  damper 
change  at  dof  13,  Ac13/0  =  2  lb-s/in  included.  The  following  figure  shows  the 
results  of  this  additional  modification. 
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Figure  5.10  Transient  Response  Beam  System  w/ Damper  Modification 

From  the  graph  produced,  it  can  been  seen  that  the  responses  are  not 
identical,  but  follow  each  other  with  relatively  little  error.  The  same  line  of 
reasoning  as  in  the  undamped  case  is  applicable  here  also.  A  smaller  time  step 
yields  more  accurate  data  but,  memory  limitations  inhibits  its  usefullness. 
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3.         Mass-Plate  System 

Now  consider  the  following  mass-plate  system: 


Akj,  =  le4  lb/in 
Ac^  =  5  lb-s/in 
Akc  =  500  lb/in 
Ac  =  5  lb-s/in 
y(t)  =  Y0(e100t -  e3001) 
Y„=l 


Figure  5.11  Plate-Mass  System  Experiencing  Base  Excitation 

The  plate  parameters  and  finite  element  modeling  of  the  above  structure  is 
identical  to  that  used  in  the  frequency  domain  synthesis  chapter  of  the  thesis. 
The  modifications  shown  are  arbitrarily  selected,  but  are  in  accordance  with 
the  changes  which  will  be  made  during  an  optimization  problem. 

a.         Spring  Modifications  Only 

The  following  figure  is  the  resultant  synthesized  and  classically 
determined  transient  response  x109(t)  which  represents  the  response  at  dof  109, 
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the  computer,  due  to  the  base  excitations  at  the  plate's  corners.  For  this 
example  none  of  the  damper  modifications  are  included. 
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Figure  5.12  Transient  Response  Mass-Plate  System 

From  the  graph  produced,  it  can  been  seen  that  the  responses  are  not  identical 
but  follow  each  other  with  relatively  little  error.  The  error  difference  between 
the  two  methods  can  be  attributed  to  the  numerical  solution  technique,  and 
the  step  size  of  the  time  vector.  As  the  step  size  decreases,  the  synthesis 
solution  approaches  the  exact  solution.  However,  memory  limitations  on  the 
computer  system  used  for  this  research,  prohibit  meaningful  use  of  a  smaller 
time  step. 
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h.  Spring  &  Damper  Modifications 

The  spring  modifications  only  example  is  now  repeated  but  with 
the  prescribed  damper  changes  at  the  base  and  between  the  computer  and  the 
plate.  The  following  figure  shows  the  results  of  these  additional 
modifications. 


0.8 


0.6 


0.4- 


E 

9 
O 

f     0 


Synthesized  Transient  Time  Response  at  dof  109 


-0.2  - 


-0.4 


-0.6 


0  0.01  0.02  0.03  0.04  0.05  0.06  0.07 

time  (sec) 


Figure  5.13  Transient  Response  Mass-Plate  System  w/ Spring-Damper 

Modifications 

Due  to  the  size  of  this  model  and  the  damper  modifications  being  made,  it 
was  not  possible  to  compare  the  synthesis  solution  with  a  classical  solution. 
The  computational  time  required  for  a  classical  solution  would  be  great. 
However,  from  comparing  Figs.  5.11  &  5.12,  the  slight  effect  of  the  additional 
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damping  can  be  seen  in  the  peak  amplitude  of  the  response,  giving  some 
credence  to  the  solution.  Also,  the  other  previously  run  damper  modification 
examples  show  that  this  technique  does  in  fact  work. 

4.         Computational  Time  Comparisons 

The  following  table  is  a  compilation  of  all  the  computation  times  that 
it  took  for  the  previously  presented  examples  to  complete  its  process. 


Synthesis 

Classical  Methods 

(sees) 

Convolution 
(sees) 

Integration 
(sees) 

At 
(sees) 

Mass- 
Spring 

No  Damper 
No  Mass 

15.62 

1.28 

.001 

Damper 
No  Mass 

23.12 

379.32 

.001 

Mass 
No  Damper 

194.05 

1.27 

.001 

Mass 
Damper 

197.89 

339.55 

.001 

Beam 

No  Damper 

107.38 

8.56 

.0003 

Damper 

183.75 

5687.88 

.0003 

Mass-Plate 

No  Damper 

27.39 

15.77 

.0008 

Damper 

99.99 

t   — »   CO 

.0008 

Table  5.1  Synthesis  vs.  Classical  Computational  Times 
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From  the  table,  it  can  be  observed  that  for  smaller  structures,  the  classical 
method  seems  more  efficient.  This  is  due  to  the  fact  that  the  advantage  of  the 
iteration  method  is  lost  when  dealing  with  small  structures.  The 
computational  time  for  the  iteration  solution,  or  in  other  words  the  synthesis 
technique,  is  driven  by  the  magnitude  of  the  time  step  and  the  number  of 
time  steps.  On  the  other  hand,  the  classical  methods'  solution  times  are 
driven  by  the  model  size.  Therefore,  depending  on  the  time  step  size,  a  very 
small  and  very  simple  model  may  take  a  very  long  time  to  solve  using  time 
domain  synthesis. 

Another  interesting  aspect  of  the  comparison  is  how  sensitive  the 
synthesis  technique  is  to  mass  changes.  Since  the  mass  change  does  not  add 
additional  dofs  to  the  model,  the  classical  method's  computational  time  is  not 
increased  by  this  modification.  On  the  other  hand,  a  mass  change  for  the 
synthesis  method  constitutes  the  use  of  the  calculated  second  derivative  of 
the  synthesized  response.  This  derivative  is  accomplished  using  a  finite 
differencing  scheme  and  inherently  adds  the  possibility  of  additional 
instability  to  the  solution. 

As  previously  alluded  to,  the  finite  differencing  could  be  the  reason 
why  the  synthesis  method  is  greatly  affected  by  mass  changes.  However, 
damper  changes,  which  also  prompt  the  use  of  finite  differencing,  did  not 
substantially  increase  computational  time  for  the  synthesis  method.  This  is  in 
contradiction  to  the  thought  that  the  increased  instability  is  caused  by  the 
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finite  differencing  scheme.  Another  factor  which  plays  into  the  convergence 
and  stability  of  the  solution  is  the  magnitude  of  the  modification.  Relative  to 
each  other  and  this  problem,  the  mass  change  may  be  of  greater  magnitude 
than  the  damper  change  and  therefore  costs  more  in  computational  time. 
The  largest  and  most  noticeable  disparity  is  in  the  time  required  to 
classically  calculate  the  dynamic  responses  following  a  damper  modification. 
The  direct  integration  method  is  very  inefficient  and  becomes  basically 
unusable  for  very  large  complex  models,  or  in  the  use  of  an  optimization 
routine. 
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VI.  OPTIMIZATION 

The  following  background  information  and  the  motivation  for  this 
chapter  can  be  attributed  in  part  to  reference  [7].  Inherent  to  the  operation  of 
military  vehicles,  whether  they  operate  on  land,  sea  or  in  the  air,  is  the 
presence  of  a  shock  and  vibration  environment.  Computer  and  electronic 
equipment  located  on  these  vehicles  can  malfunction  or  be  damaged  as  a 
result  of  the  severe  effects  of  this  environment.  Therefore,  the  need  arises  for 
a  properly  designed  isolation  system,  which  will  absorb  the  shock  and 
vibrations  and  therefore  protect  the  equipment.  This  is  an  especially 
important  concept  in  today's  military  with  the  relatively  recent  practice  of 
using  Commercial-Off-The-Shelf  (COTS)  equipment  vice  the  traditionally 
expensive  shock  and  vibration  hardened  military  equipment. 

In  general,  it  is  relatively  straightforward  to  design  an  isolation  system 
to  protect  against  shock  only  or  vibration  only.  However,  these  two  designs 
are  in  direct  conflict  with  each  other  [Ref.  1].  Therefore,  in  the  presence  of 
both  types  of  excitations  and  various  other  constraints,  the  optimal  design  of 
an  isolation  system  becomes  more  complicated.  The  need  then  arises  for 
some  sort  of  optimization  scheme  which  will  minimize  the  response  of  the 
component  to  be  protected,  while  simultaneously  satisfying  the  constraints 
which  have  been  placed  on  the  system. 
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The  complexity  of  this  problem  does  not  necessarily  end  here.  If  the 
structure  to  be  isolated  and  the  isolation  system  to  be  designed  is  large  and 
complex,  the  analysis  process  to  calculate  the  desired  responses  can  be  very 
computationally  demanding.  Furthermore,  if  there  are  multiple  design 
variables  which  may  be  altered  in  order  to  achieve  the  desired  result,  the 
search  patterns  for  the  optimal  solution  are  more  complex,  which  leads  to 
more  iterations  of  the  optimization  process.  The  combination  of  these  two 
factors  can  lead  to  the  conclusion  that  the  optimal  design  of  a  large  complex 
isolation  system  is  impractical. 

With  the  previously  presented  arguments  in  mind,  the  use  of 
computationally  efficient  synthesis  techniques  for  the  calculation  of  the 
system's  response  is  the  obvious  answer.  Since  the  structure  is  large  and 
complex,  analysis  time  could  be  relatively  long.  The  larger  the  system,  the 
greater  the  number  of  possible  design  variables.  This  means  that  the  search 
for  the  optimal  solution  is  even  more  complex,  and  will  require  more 
iterations.  This  portion  of  the  thesis  will  demonstrate  the  use  of  the 
previously  presented  structural  synthesis  techniques  in  the  optimal  design  of 
a  shock  and  vibration  isolation  system.  A  general  computer  algorithm  will  be 
formulated  and  used  in  order  to  illustrate  the  use  of  structural  synthesis  in 
optimization  design. 


82 


A.        GENERAL  OPTIMIZATION  PROBLEM  FORMULATION 

In  an  optimization  problem,  the  quantity  to  be  maximized  or 
minimized  is  termed  the  objective  function.  The  quantities  which  may  be 
altered  in  order  to  find  this  optimal  value  are  termed  the  design  variables. 
The  constraint  functions  are  also  a  function  of  the  design  variables,  and  three 
different  types  of  constraints  may  exist  in  an  optimization  problem.  The  first 
type  is  an  inequality  constraint,  where  the  value  calculated  from  the 
constraint  function  is  less  than  or  equal  to  some  prescribed  limit.  The  second 
type  is  an  equality  constraint.  Here  the  value  calculated  from  the  constraint 
function  must  be  equal  to  some  prescribed  value.  The  third  constraint  type  is 
side  constraints.  This  is  when  upper  and  lower  limits  are  placed  on  the  design 
variables.  In  a  constrained  optimization  problem,  the  objective  function  is 
maximized  or  minimized,  while  ensuring  that  the  user  defined  values  of  the 
constraint  functions  are  not  violated. 

For  the  design  of  an  isolation  system,  the  minimization  of  the 
equipment's  response  over  some  frequency  range,  due  to  vibration,  could  be 
the  main  concern.  At  the  same  time,  the  entire  structure's  response  due  to  a 
shock  input,  and  the  static  displacement  or  'sag7  of  the  system  could  be 
constrained.  The  physical  parameters  which  can  be  modified  to  achieve  this 
goal  could  be  the  stiffness  and  damping  values  of  the  isolators.  The  design  of 
this  isolation  system  can  be  converted  to  the  following  optimization  problem: 
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Minimize  (Objective  Function): 

Computer  response  due  to  random  vibration 

Subject  to  (Constraints): 

Maximum  static  displacements  <  limit 
Maximum  dynamic  isolator  stretch  due  to  shock  <  limit 
Maximum  computer  acceleration  due  to  shock  <  limit 
Min  and  max  changes  in  the  isolator  stiffness  values  <  limit 
Min  and  max  changes  in  the  isolator  damping  values  <  limit 

Design  Variables: 

Changes  in  the  isolator  stiffness  values 
Changes  in  the  isolator  damping  values 

1.         Calculation  of  Objective  Function 

Since  the  objective  function  is  based  on  the  response  due  to  random 
vibrations,  the  input  base  excitation  is  in  the  form  of  a  power  spectral  density 
(PSD)  function.  The  response  is  therefore  calculated  as  the  mean  square  value 
in  terms  of  the  system  response  function  (FRF),  and  the  PSD  of  the  input.  The 
equation  for  the  mean  square  value  of  the  response  is  given  as  [Ref.  10]: 

x2  =  £|H(/)fs(/)d/  (6.D 

From  equation  (6.1),  it  can  be  seen  that  the  most  accurate  means  of 
minimizing  this  equation  would  be  if  the  input  PSD,  S(/),  is  known.  Since  the 
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goal  is  to  formulate  a  general  optimization  program,  it  is  not  inconceivable 
and  is  appropriate  to  assume  S(/)  to  be  constant  over  some  frequency,  fl<f< 
f2.  Therefore,  minimizing  the  mean  square  response  over  the  appropriate 
frequency  range  is  equivalent  to  minimizing: 

F  =  J/2|H(/)|2d/  (6.2a) 

It  can  be  seen  that  for  every  iteration  of  the  optimization  process,  the 
calculations  involved  in  equation  (6.2a)  must  be  performed.  If  the  FRFs  of  the 
baseline  structure  is  provided  or  calculated  prior  to  the  iterations  begin,  the 
frequency  domain  structural  synthesis  technique  may  be  used  to  calculate  all 
subsequent  FRFs  and  the  objective  function  is  properly  written  as: 

F  =  J^|H'(/)|2d/  (6.2b) 

2.         Calculation  of  Constraints 

The  first  constraints  mentioned  are  the  maximum  static  displacements 
of  the  structure.  Since  these  values  must  also  be  recalculated  every  iteration, 
an  efficient  means  of  doing  so  needed.  Various  reasons  make  this  problem  a 
prime  candidate  for  the  use  of  static  displacement  synthesis: 

•  the  static  displacement  at  select  locations  vice  the  entire  structure  is 
desired,  therefore  calling  for  some  type  of  reduction  or 
condensation  technique. 

•  the  elemental  and  system  matrices  of  the  structure  may  not  be 
available. 
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•  the  modal  solution  of  the  eigenvalue  problem  has  already  been 
accomplished. 

•  the  isolator  stiffness,  which  are  design  variables,  act  in  vertical 
translation  direction  only  (isolator  damping  has  no  effect  on  static 
problem). 

Another  constraint  is  the  maximum  dynamic  isolator  stretch  due  to 
the  shock  input.  In  order  to  calculate  this  constraint,  the  transient  response  of 
the  structure  at  the  dofs  where  the  isolators  are  connected  must  be  obtained. 
The  displacement  of  the  base  input  as  a  function  of  time  must  also  be  known. 
Since  the  change  in  the  isolator  damper  represents  nonproportional 
damping,  the  direct  integration  method  would  be  used  to  solve  this  problem. 
This  calculation  would  be  terrible,  in  reference  to  computing  time  for  a  single 
iteration  of  the  optimization  process.  Despite  other  reasons,  this  alone  is 
enough  motivation  for  the  implementation  of  the  time  domain  synthesis 
technique  for  the  calculation  of  the  dynamic  responses.  The  last  constraint, 
the  computer  acceleration,  is  obtained  by  simply  applying  a  finite  differencing 
scheme  to  the  computer's  displacement  time  history.  All  the  above 
mentioned  constraints  are  known  as  inequality  constraints,  and  will  satisfy 
the  general  normalized  form  gj(X)  <  0  [Ref.  11]. 
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3.         Design  Variables 

There  are  limits  on  the  amount  in  which  the  isolator  stiffness  and 
damping  can  be  changed.  For  example,  the  lower  bound  for  the  change  in  an 
isolator  stiffness  should  not  be  equal  to  the  negative  of  the  original  value  of 
the  isolator.  If  it  is,  the  optimization  program  might  use  this  value  and  this 
constitutes  the  total  removal  of  the  isolator  and  changing  the  basic  structure 
of  the  system.  The  upper  bound  for  the  change  is  limited  by  what  is  physically 
realistic.  These  type  of  constraints  on  the  design  variables  are  known  as  side 
constraints  [Ref.  11]. 

The  relative  value  of  design  variables  to  one  another  is  very  important 
when  considering  the  efficiency  and  reliability  of  the  numerical  optimization 
process  [Ref.  11].  The  function  space  of  the  problem,  which  is  defined  by  the 
design  variables,  should  be  as  symmetric  as  possible.  When  looking  at  the 
design  variables,  isolator  stiffness  and  isolator  damping,  it  can  be  seen  that 
effective  changes  to  isolator  stiffness  are  on  the  magnitude  of  thousands, 
while  the  isolator  damper  changes  are  on  the  magnitude  of  ones.  This  large 
disparity  creates  a  very  nonsymmetric  function  space.  In  order  to  alleviate 
these  problems,  the  design  variables  are  scaled  so  that  their  magnitudes  are 
comparable.  A  symmetric  function  domain  allows  for  the  faster 
determination  of  the  optimization  by  decreasing  the  number  of  iterations 
necessary  to  find  the  optimal  solution  [Ref.  11]. 
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B.         OPTIMIZATION  COMPUTER  CODE 

The  optimization  tool  used  to  solve  this  problem  is  the  MATLAB 
optimization  toolbox  l.Od.  More  specifically,  the  function  constr.m,  which  is 
designed  for  optimization  of  a  constrained,  nonlinear,  multivariable 
problem,  is  used.   This  function  uses  the  Sequential  Quadratic  Programming 
(SQP)  method  in  order  to  solve  this  problem  [Ref.  12].  The  SQP  method  is 
used  in  order  to  determine  search  directions  for  the  design  variables  at  every 
iteration  [Refs.  11  &  12].  One  of  the  limitations  of  this,  and  other  traditional 
nonlinear  optimization  codes,  is  that  they  may  only  give  local  solutions. 
Also,  when  the  problem  is  infeasible,  constr.m  attempts  to  minimize  the 
maximum  constraint  value.  What  this  can  lead  to  is,  the  problem  iterating  to 
the  user  defined  iteration  limit,  the  search  pattern  extending  outside  the 
bounds  of  the  design  variables,  and  no  optimum  solution  ever  being  found. 
Therefore,  the  user  must  have  good  understanding  of  the  problem  and  the 
constraints  imposed  upon  it. 

The  goal  in  the  formulation  of  the  optimization  computer  programs  is 
to  provide  a  means  of  determining  the  optimal  changes  in  isolator  values,  in 
order  to  minimize  some  dynamic  response,  while  simultaneously  satisfying 
all  prescribed  constraints.  These  programs  first  allow  for  the  loading  of  a 
general  finite  element  model  of  the  structure.  This  baseline  structure  is  given 
in  terms  of  its  FRFs  and  IRFs,  or  in  terms  of  its  system  matrices  from  which 
the  baseline  FRFs  and  Ifs  must  be  calculated.  Modifications  to  the  baseline 
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structure  which  will  act  as  design  variables  are  then  designated.  User  defined 
initial  values  of  design  variables  and  constraint  values  are  designated.  During 
the  optimization  process,  all  structure  responses  will  be  calculated  using 
structural  synthesis  techniques.  With  minor  coding  adjustments,  the  ability 
to  interchange  objective  function  and  constraints  exists.  For  example  the 
optimization  program  could  be  the  minimization  of  the  computer's 
acceleration  due  to  shock,  while  the  response  to  a  random  vibration  input  is  a 
constraint.  Appendix  D  shows  the  computer  codes  used  to  perform  this 
process.  The  following  figure  is  a  diagram  of  the  flowpath  of  the  optimization 
program. 
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Figure  6.1  Flowpath  of  the  Optimization  Program 
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C         ILLUSTRATIVE  EXAMPLE 

The  purpose  of  this  section  is  to  illustrate  the  use  of  the  structural 
synthesis  techniques  in  the  optimization  of  a  shock  and  vibration  isolation 
system.  The  figure  below  shows  the  mass-plate  structural  system  which  was 
used  in  the  previous  examples  and  is  also  used  here. 


y(t)f 


Figure  6.2  Plate-Mass  System  Experiencing  Base  Excitation 


1.         Problem  Formulation 

The  following  is  the  formal  problem  statement  of  the  optimization 
problem: 

r!60,  .2 


Minimize: 


F(Akb,Akc,Acb/Acc)  =  Jo    |H*(/)|d/ 
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Subject  to: 


gi  = 


g2  = 


g3  = 


g4  = 


g5  = 


z„  stat 

-*= 1<0 


max  zp_stat 

zc_stat 
max  z„  stat 


-1<0 


k„  stretch 

-1<0 


max  kp_stretch 

kc_stretch 
max  k„  stretch 


1<0 


zc_accel 
max  zc_accel 


-1<0 


where: 


^000  <Akb<  5000 
-4000  <Ak<  5000 


0  <  Acb  <  5 
0  <  Acc  <  5 


H*(/)  =  frequency  response  function  of  the  computer 

Zp_stat  =  plate  static  displacement  [in] 

zc_stat  =  relative  static  displacement  between  computer  and  plate  [in] 

maxZp_stat  =  maximum  allowable  plate  static  displacement  [in] 

maxzc_stat  =  maximum  allowable  static  displacement  between 

computer  and  plate  [in] 
kp_stretch  =  isolator  stretch  between  base  and  plate  due  to  shock  [in] 
k^stretch  =  isolator  stretch  between  plate  and  computer  due  to  shock 

[in] 
max  kp_stretch  =  maximum  allowable  isolator  stretch  between  base  and 

plate  [in] 
max  kc_stretch  =  max  allowable  isolator  stretch  between  plate  and 

computer  [in] 
z^accel  =  absolute  acceleration  of  computer  due  to  shock  [in/s2] 
max  zc_accel  =  max  allowable  absolute  acceleration  of  computer  [in/s2] 
Akb/  A  kc,  Acb/  Acc  =  the  design  variables. 
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The  following  is  information  that  must  be  provided  to  the 
optimization  program  for  its  execution,  and  for  the  understanding  of  the 
user.  The  following  table  lists  the  initial,  lower  bound,  and  upper  bound 
values  for  the  design  variables. 


Ak> 


Ak 


Ac, 


Ac, 


initial  value 
lower  bound 
upper  bound 


1000 

-4000 

5000 


1000 
4000 
5000 


Table  6.1  Optimization  Inputs 

From  the  table  it  can  be  noticed  that  an  initial  value  for  the  design  variables  is 
used.  Since  the  optimization  routine  locates  local  minimum,  the  start 
position  of  the  search  is  very  important.  The  baseline  values  for  the  isolator 
elements  which  are  to  be  modified  are: 


initial  kb  =  10000  lb /in 
initial  kc  =  5000  lb/in 
initial  cb  =0  lb-s/in 
initial  cc  =0  lb-s/in 


The  constraint  values  for  the  inequality  constraint  are: 

maxz    stat  =  -0.08  in 

p— 

maxz^stat  =  -0.03  in 
max  kp_stretch  =  1.0  in 
max  kc_stretch  =  1.0  in 
max  zc_accel  =  30g's 
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2.         Optimization  Results 

The  previously  described  problem  is  executed  using  the  isomain.m  and 
isosub.m  programs.  During  the  execution  of  the  optimization  program, 
output  information  is  provided  through  the  MATLAB  Command  Window 
(Table  6.2). 


f-COUNT 

FUNCTION 

MAXfg) 

STEP 

Procedures 

5 

36012.1 

1.40603 

1 

10 

15494 

0.0364448 

1 

15 

6969.56 

-0.0262839 

1 

22 

6691.27 

-0.0370827 

0.25 

29 

6372.67 

-0.0682749 

0.25 

34 

6372.59 

-0.0693421 

1 

39 

6370.5 

-0.0631356 

1 

mod  Hess(2) 

44 

6370.49 

-0.0687934 

1 

mod  Hess 

49 

6370.49 

-0.0689084 

1 

mod  Hess(2) 

54 

6370.48 

-0.067107 

1 

59 

6370.47 

-0.0527099 

1 

mod  Hess 

64 

6370.45 

-0.0667607 

1 

69 

6370.44 

-0.0687255 

1 

mod  Hess(2) 

74 

6368.61 

-0.0685718 

1 

mod  Hess(2) 

79 

6368.61 

-0.0688113 

1 

mod  Hess 

80 

6368.61 

-0.0667838 

1 

mod  Hess 

Table  6.2  Optimization  l.Od  Generated  Output 

Once  the  optimization  program  terminates  successfully,  post 
processing  occurs  and  the  following  figures  and  summary  information  is 
provided.  The  first  figure  is  a  graph  of  the  value  of  the  design  variables  at 
each  iteration. 
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Figure  6.3  Values  of  the  Design  Variables 

The  next  figure  is  a  graph  of  the  absolute  value  of  the  isolators  at  each 
iteration  of  the  optimization  process. 
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Change  in  Isolation  Constants  vs.  Iterations 
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Change  in  Isolation  Constants  vs.  Iterations 
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Figure  6.4  Absolute  Values  of  the  Isolator  Elements 

The  next  figure  is  a  graph  of  the  FRF  of  the  computer  before  the  optimization 
begin,  and  at  its  termination. 
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Optimized  Frequency  Response  Function  [H*] 


60  80  100  120  140  160 

Frequency  (Hz) 


Figure  6.5  Computer  FRF  Before  and  After  Optimization 

The  effect  of  the  changes  to  the  structure,  particularly  the  damping  additions, 
can  definitely  be  seen  in  this  figure.  The  objective  function  is  defined  as  the 
area  under  this  curve  after  it  is  squared.  Therefore,  by  lowering  the  peak 
values  of  the  response,  the  objective  function  is  ultimately  lowered.  The  next 
figure  is  a  plot  of  how  the  value  for  the  objective  function  changed  as  a 
function  of  the  optimization  iterations. 
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Figure  6.6  Values  of  the  Objective  Function 

From  Fig.  6.6,  the  dramatic  effects  of  the  optimization  can  really  be  seen.  The 
initial  value  of  the  objective  function  was  36012.1  and  after  the  isolator 
modifications  the  objective  function  equals  6368.31.  The  following 
information  is  a  summary  of  the  optimization  process  and  is  displayed  at  the 
MATLAB  Command  Window: 


Optimization  Terminated  Successfully 

Active  Constraints: 

ans  =  3 

The  recommended  change  in  base  isolator  stiffness  (lb /in)  is 
final_del  kb  =  2.3271e+03  lb/in 


— > 
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The  recommended  change  in  computer  isolator  stiffness  (lb/ in)  is  — > 
final_del_kc  =  4.0849e+03  lb/in 


The  recommended  change  in  base  isolator  damping  (lb/in/s)  is  — > 
final_del_cb  =  2.7678  lb-s/in 


The  recommended  change  in  computer  isolator  damping  (lb/in/s)  is  — > 
final_del_cc  =  4.6266  lb-s/in 

Constraints: 

gj  =  -0.1814  g2  =  -0.6334  g3  =  -0.5237  g4  =  -0.6981  g5  =  -0.0668 

elapsed_time  =  9.5550e+03  sees  =  2hrs  39min 


It  can  be  seen  from  the  presented  information  that  it  took  the 
optimization  program  80  iterations  to  recommend  the  presented  changes  to 
the  structure  in  order  to  optimize  its  response  against  random  vibrations,  and 
meet  the  established  constraint  criteria.  The  most  notable  portion  of  this 
example  is  that  it  accomplished  this  task  on  a  109  dof  structure,  which 
underwent  non-proportional  damping  changes  at  every  iteration,  in  only  2 
hours  and  39  minutes.  In  other  words,  this  program  performed  80  transient 
response  analyses,  calculated  the  static  displacement  80  times,  and  calculated 
the  frequency  response  function  and  the  integral  of  that  curve  80  times,  for  a 
109  dof  structure. 
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VII.  CONCLUSIONS  /RECOMMENDATIONS 

The  main  objective  of  this  thesis  was  to  investigate  the  use  of 
structural  synthesis  techniques  as  rapid  reanalysis  procedures  in  the  optimal 
design  of  large  shock  and  vibration  isolation  systems.  In  order  to  use  these 
synthesis  techniques  liberally,  an  investigation  into  the  accuracy  of  the 
methods,  and  the  computer  algorithms  formulated,  had  to  be  accomplished. 
It  is  important  to  mention  at  this  point  that  the  time  and  frequency  structural 
synthesis  methods  used  in  this  study,  primarily  covers  base  excitations.  This 
is  a  new  extension  of  the  more  familiar  synthesis  formulations.  The  static 
displacement  synthesis,  is  however  a  completely  new  development.  These 
formulations  had  to  be  general  enough  to  handle  different  structures  and  to 
handle  various  types  of  modifications  to  the  structure.  The  synthesis 
techniques  were  then  coupled  with  an  optimization  routine,  and  together 
they  were  used  in  finding  the  optimal  changes  to  a  given  isolation  system. 

A.        CONCLUSIONS 

From  the  various  analyses  conducted,  the  following  conclusions  can  be 
drawn: 

1.  The  frequency  domain  structural  synthesis  technique  is  a  very 
efficient  and  accurate  method  of  calculating  the  new  FRFs  of  a 
structure  following  indirect  modifications.  It  is  nonrestrictive  and 
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arbitrary  as  to  where  and  what  type  of  modifications  can  be  made  to 
the  structure. 

2.  The  static  displacement  structural  synthesis  is  also  a  very  efficient 
and  accurate  method  of  calculating  the  new  static  displacement  of  a 
structure  following  modifications.  The  accuracy  of  this  method  does 
depend  on  the  use  of  a  finite  differencing  scheme  which  in  turn  is 
dependent  on  the  variable  step  size.  It  is  also  found  that  this  method 
is  restrictive  in  its  use.  Certain  criteria,  which  are  identified  in 
Chapter  IV  of  this  work,  must  be  met  in  order  that  this  synthesis 
technique  can  be  applied  properly.  For  the  isolation  problems 
presented  in  this  work,  the  required  criteria  were  met. 

3.  The  time  domain  structural  synthesis  technique  is  an  efficient  and 
accurate  method  of  calculating  the  new  transient  response  of  a 
structure  following  indirect  modifications.  It  is  nonrestrictive  and 
arbitrary  as  to  where  and  what  type  of  modifications  can  be  made  to 
the  structure.  However,  since  numerical  techniques  are  used  to 
solve  for  the  responses,  some  error  of  approximation  is  entered  into 
the  solution.  As  the  time  step  is  decreased,  the  numerical  solution 
approaches  the  exact  solution.  However,  current  solution 
algorithms  are  restricted  by  memory  requirements. 

4.  In  finding  the  optimal  design  of  a  relatively  large  shock  and 
vibration  isolation  system,  the  use  of  the  synthesis  techniques 
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proved  to  be  an  invaluable  asset  in  the  optimization  routine.  With 
the  use  of  the  synthesis  techniques,  in  2  hours  and  39  minutes,  it 
was  possible  to  perform  80  FRF,  static  displacement,  and  transient 
analyses  on  a  109  dof  system,  which  was  non-proportionally 
damped.  To  get  an  idea  of  how  efficient  this  is,  recall  from  the  time 
domain  synthesis  analysis  (Chapter  V),  that  it  took  approximately  1 
and  1/2  hours  to  calculate  one  transient  analysis  of  a  22  dof  beam 
structure. 

B.         RECOMMENDATIONS 

Throughout  this  study,  some  weaknesses  of  the  techniques,  and 
opportunities  to  enhance  the  techniques  were  touched  upon. 
Recommendations  to  address  these  weaknesses,  and  take  advantage  of  these 
opportunities  include: 

•  All  modal  analyses  were  done  to  include  all  modes  of  the  structure. 
A  modal  truncation  criteria  should  be  established,  and  studies 
should  be  done  to  ensure  that  these  criteria  hold  up  for  synthesis 
also.  This  way  it  would  be  possible  to  take  advantage  of  time  saving 
aspects  which  are  inherent  to  modal  superposition  techniques. 

•  Use  actual  FRF  and  IRF  data  from  a  real  structure,  and  predict 
structural  responses  due  to  modifications. 
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•  Build  CAD  models  to  validate  damping  modifications  to  large 
structures. 

•  Investigate  to  obtain  a  general  rule  for  choosing  the  optimum  AQ,  to 
be  used  in  finite  difference  scheme  to  calculate  static  response. 
Another  option  would  be  to  investigate  the  use  of  the  closed  form 

2nd  derivative  of  H*   to  obtain ^5 — . 

da2 

•  Implement  relaxation  schemes  for  the  iterative  solution  in  time 
domain  synthesis.  Also  a  study  should  be  done  in  order  to 
determine  what  the  proper  time  step  size  should  be  to  ensure 
convergence  and  stability. 
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APPENDIX  A.  FREQUENCY  DOMAIN  SYNTHESIS  COMPUTER 

CODES 

The  following  is  a  list  and  a  brief  description  of  the  main  MATLAB 
computer  codes  that  were  written  in  order  to  perform  the  frequency  domain 
synthesis  calculations. 

•  fsynstr.m  -  performs  the  frequency  domain  synthesis  on  a  structure  which 
experiences  an  external  force  excitation  directly  on  the  structure,  and 
compares  this  solution  to  the  FRF  calculated  using  classical  methods. 

•  fsynbase.m  -  performs  the  frequency  domain  synthesis  on  a  structure 
which  experiences  a  base  excitation,  and  compares  this  solution  to  the  FRF 
calculated  using  classical  methods. 

•  fsynmain.m  -  the  main  program  which  calls  the  modularized  frequency 
synthesis  programs  and  does  the  post  processing  of  the  results. 

•  change.m  -  loads  the  baseline  structure  to  be  modified,  allows  for 
modifications  to  the  baseline  structure,  and  calls  the  function  delta.m 

•  delta.m  -  forms  the  modification  matrices  which  will  be  used  to  form 
impedance  matrices  or  determine  coupling  forces. 

•  modal.m  -  solves  for  the  structure's  natural  frequencies,  mode  shapes  and 
other  modal  matrices  and  vectors. 

•  frfmodal.m  -  calculates  the  structure's  FRFs. 
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•  frfsynth.m  -  performs  the  synthesis  on  the  baseline  structure  and  returns 
the  synthesized  FRFs. 

•  frf  sift.m  -  selects  the  FRF  that  the  user  wants  to  perform  post  processing 
on. 

The  full  codes  are  contained  on  subsequent  pages. 
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fsynstr.m 


%  This  program  is  designed  to  calculate  FRFs  by  inverting  the  impedance  matrix. 

%  Changing  the  mass,  stiffness,  and  damping  matrices,  new  FRFs  will  be  calculated  by 

%  reformulating  the  matrices,  and  using  the  synthesis  method  in  the  frequency  domain. 

% 

clear 

%  Formation  of  the  elemental  matrices 

% 

nel=4; 

nn=nel+l; 

dof=nel+l; 

ndk=l; 
knel(l)=0; 
for  spring=l:ndk; 
kcon(spring)= 100; 
knel(spring+ 1  )=4; 

for  ele=knel(spring)+ 1  :knel(spring+ 1 ) 

kcon_ele(ele)=kcon(spring); 
end 
end 

ndm=l; 

mnel(l)=0; 

formass=l:ndm; 

m(mass)=100/386.4; 
mnel(mass+l)=5; 

for  ele=mnel(mass)-i-l  :mnel(mass+l) 

mele(ele)=m(mass); 
end 
end 

%  Formulation  of  global  stiffness,  damping  and  mass  matrices  from  elemental  matrices 
% 

Kglo=zeros(dof,dof); 
Cglo=zeros(dof,dof); 
Mglo=zeros(dof,dof); 
forn=l:nel 
kele=kcon_ele(n)*[  1-1 

-i  l]; 

rcl=n; 

rc2=n+l; 

Kglo(rcl:rc2^-cl:rc2)=Kglo(rcl:rc2jcl:rc2)  +  kele; 
end 
Mglo=diag(mele); 
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alpha=0;    beta=.0005; 
Cglo=alpha*Mglo  +  beta*Kglo; 

bc=2; 
ifbc=l 

K=Kglo(2:dof,2:dof); 

K((dof-l),:)=[]; 

K(:,(dof-l))=[]; 

M=Mglo(2:dof,2:dof); 

M((dof-l),:)=[]; 

M(:,(dof-l))=[]; 

C=Cglo(2:dof,2:dof); 

C((dof-l),:)=[]; 

C(:,(dof-l))=[]; 

ndof=dof-2; 
elseif  bc=2 

K=Kglo(2:dof,2:dof); 

M=Mglo(2:dof,2:dof); 

C=Cglo(2:dof,2:dof); 

ndof=dof-l; 
elseif  bc=3 

K=Kglo(2:dof,2:dof); 

K((dof-l),:)=[]; 

K(:,(dof-l))=[]; 

M=Mglo(2:dof,2:dof); 

M((dof-l),:)=[]; 

M(:,(dof-l))=[]; 

C=Cglo(2:dof,2:dof); 

C((dof-l),:)=[]; 

C(:,(dof-l))=[]; 

ndof=dof-2; 
else 

K=Kglo;       M=Mglo;  C=Cglo; 

ndof=dof; 
end 

%  This  portion  of  the  program  is  designed  to  accept  changes  to  the  original 
%  structure  and  recalculate  the  FRF  using  classical  method  of  reforming  [M], 
%  [K],  and  [C]  matrices  using  the  inverse  of  the  impedence  matrix.  Also 
%  recalculates  the  FRF  using  structural  synthesis  in  the  frequency 
%  domain. 
% 

change_model=  1 ; 
while  change_model=l; 

K_c=K;     M_c=M;      C_c=C;  %  initializes  change  matrices  to 

%  original  matrices 

changek=2;    changec=2;     changem=2;     %  initializes  logic  variables 

%  to  'no'  status 
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chk=0;    chc=0;     chm=0;     %  initializes  counter  for  #  of  changes 

%  to  matrices 

kdof  1=[];      kdof2=D;         %  initalizes  the  matrices  which  keep 
cdof  1=[];      cdof2=Q;         %  track  of  the  dofs  where  changes  occur 
mdof=D; 

cset=0;  %  initializes  cset  matrix  to  zero  to  avoid  null  index  error  in 

%  handling  the  spring  to  ground  if  no  changes  are  made  to  the 
%  structure. 

changek=inputCWould  you  like  to  change  your  stiffness  matrix?  l=yes  2=no  ')» 
while  changek=l 

chk=chk+l;     %  counter  for  determining  #  changes  to  stiffness  matrix. 

lk(chk)=input('input  value  of  added  spring  (lbs/in)  '); 

kdofl(chk)=input('input  constrained  dof  where  1st  end  of  added  spring  is  applied  '); 

kdof2(chk)=input('input  dof  where  2nd  end  of  added  spring  is  applied  or  0  for  ground 

'); 

%  comparison  of  change  dof  with  the  c-set  matrix  and  formation  of 

%  new  c-set  matrix. 

% 

if  kdof  1  (chk)~=cset 

cset=[cset  kdofl(chk)]; 
end 
if  kdof2(chk)~=cset 

cset=[cset  kdof2(chk)]; 
end 

%  Classical  method  of  reformulatting  the  [K]  matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk))  +  lk(chk); 

ifkdof2(chk)~=0 

K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)4cdof2(chk))  -  lk(chk); 

K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)Jcdofl(chk))  -  lk(chk); 

K_c(kdof2(chk)  Jcdof2(chk))=K_c(kdof2(chk)^dof2(chk))  +  lk(chk); 
end 

changek=input('Are  there  any  other  changes  to  the  stiffness  matrix?  l=yes  2=no  '); 

end 

changec=inputCWould  you  like  to  change  your  damping  matrix?  l=yes  2=no  '); 
while  changec=l 

chc=chc+l ;     %  counter  for  determining  #  changes  to  damping  matrix. 

flag_c=l;  %  mapping  matrix  flag 

lc(chc)=input('input  value  of  added  damper  (lbs-sec/in)  '); 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is  applied  '); 

cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied  or  0  for 
ground  '); 
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%  comparison  of  change  dof  with  the  c-set  matrix  and  formation  of 

%  new  c-set  matrix. 

% 

if  cdof  1  (chc)~=cset 

cset=[cset  cdofl(chc)]; 
end 
if  cdof2(chc)~=cset 

cset=[cset  cdof2(chc)]; 
end 


%  Classical  method  of  reformulatting  the  [C]  matrix. 

% 

C_c(cdofl(chc),cdofl(chc))=C_c(cdofl(chc),cdofl(chc))  +  lc(chc); 

ifcdof2(chc)~=0 

C_c(cdofl(chc),cdof2(chc))=C_c(cdofl(chc),cdof2(chc))  -  lc(chc); 

C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc))  -  lc(chc); 

C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc))  +  lc(chc); 
end 

changec=input('Are  there  any  other  changes  to  the  damping  matrix?  l=yes  2=no  '); 
end 

changem=inputCWould  you  like  to  change  your  mass  matrix?  l=yes  2=no  '); 

while  changem=l 
chm=chm+l ;         %  counter  for  determining  #  changes  to  mass  matrix. 
lm(chm)=input('input  value  of  added  mass  (lbs-secA2/in)  ')', 
mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 

%  comparison  of  change  dof  with  the  c-set  matrix  and  formation  of 

%  new  c-set  matrix. 

% 

if  mdof (chm)~=cset 

cset=[cset  mdof(chm)]; 
end 

%  Classical  method  of  reformulatting  the  [M]  matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm))  +  lm(chm); 

changem=dnput('Are  there  any  other  changes  to  the  mass  matrix?  l=yes  2=no  '); 
end 

%  Handling  of  spring  added  to  ground 

% 

ground  =  find(cset=0); 

cset(ground)  =  Q; 


110 


%  Formulation  of  iset  matrix  which  is  the  set  of  dofs  other  than 

%  those  in  the  cset  matrix  where  response  information  is  desired. 

% 

iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in?  '); 

%  Formation  of  eset  matrix  which  is  i  U  c. 

% 

eset=[iset  cset]; 

%  Choosing  the  FRF  interested  in 
% 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

inp_dof=inputCWhat  input  dof  are  you  interested  in?  '); 

frf_index=find(eset=resp_dof  I  eset=inp_dof); 

if  resp_dof =inp_dof 
frf_index=[frf_index  frf_index]; 

end 
%  Structural  Synthesis  method  of  calculating  FRF 
% 

%  Calculation  of  unchanged  FRF 

% 

freq=.01:.01:10; 

omega=2*pi*freq; 

j=length(omega); 

for  w=l:j 

Z=K+i*omega(w)*C-omega(w)A2*M; 

H=inv(Z); 

%  Rearrangement  of  original  [H]  to  form  [Hee],  [Hec],  [Hcc],  [Hce] 

%  Also  provides  for  if  user  makes  no  changes  to  the  model 

% 

Hee=H(eset,eset); 

if  cset=[] 

Hec=0;  Hcc=0;  Hce=0; 

else 

Hec=H(eset,cset); 

Hcc=H(cset,cset); 

Hce=H(cset,eset); 
end 

%  Formation  of  [del_K],  [del_C],  [del_M],  &  [del_Z] 

% 

del_K=zeros(length(cset));  %  intializes  the  change  matrices 

del_C=zeros(length(cset));  %  to  zero  and  sets  their  sizes 

del_M=zeros(length(cset));  %  as  [cset  x  cset]  matrix 

ifchk~=0 
for  a=l:chk 

indexkl=find(kdofl(a)==cset);    %  finds  where  in  cset  matrix 
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indexk2=find(kdof2(a)=cset);    %  the  change  in  k  occurs 

del_K(indexkl,indexkl)=del_K(indexkl4ndexkl)  +  lk(a); 
del_K(indexkl,indexk2)=del_K(indexkl,indexk2)  -  lk(a); 
del_K(indexk2,indexkl)=del_K(indexk2,indexkl)  -  lk(a); 
del_K(indexk2,indexk2)=del_K(indexk2,indexk2)  +  lk(a); 
end 
end 

if  chc~=0 
fora=l:chc 
indexc  1  =find(cdof  1  (a)==cset);    %  finds  where  in  cset  matrix 
indexc2=find(cdof2(a)=cset);    %  the  change  in  k  occurs 

del_C(indexcl, indexc  l)=del_C(indexcl4ndexcl)  +  lc(a); 
del_C(indexcl,indexc2)=del_C(indexcl,indexc2)  -  lc(a); 
del_C(indexc2,indexcl)=del_C(indexc2,indexcl)  -  lc(a); 
del_C(indexc2,indexc2)=del_C(indexc2,indexc2)  +  lc(a); 
end 
end 

if  chm~=0 
fora=l:chm 
indexm=find(mdof(a)==cset);     %  finds  where  in  cset  matrix 

%  the  change  in  k  occurs 

del_M(indexm,indexm)=del_M(indexm,indexm)  +  lm(a); 
end 
end 

ifcset==[] 

del_Z=0; 
else 

del_Z=del_K+i*omega(w)*del_C-omega(w)A2*del_M; 
end 

H_star=Hee-Hec*del_Z*inv(eye(size(Hcc,l))  +  Hcc*del_Z)*Hce; 

%  Capturing  the  FRF  interested  in  at  each  frequency  swept  thru 
% 
hstar(w)=H_star(frf_index(  1  ),frf_index(2)); 


%  Classical  method  of  calculating  the  FRF  by  reassembling  [Z]  and  taking  the 
%  inverse. 

% 

Z_c=K_c+i*omega(w)*C_c-omega(w)A2*M_c; 

H_c=inv(Z_c); 

H_c=H_c(eset,eset); 
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%  Capturing  the  FRF  interested  in  at  each  frequency  swept  thru 

% 

h_c(w)=H_c(ftf_index(  1  ),frf_index(2)); 


end 


%  Plotting  of  FRF  using  Classical  &  Frequency  Synthesis  Methods 
% 

semilogy(freq,abs(h_c),'r— ',fTeq,abs(hstar),'b') 

title([ frequency  Response  Function  (H,,int2str(resp_dof),int2str(inp_dof),1)'  ]) 

ylabel(*FRF  Amplitude  (in/lb)') 

xlabelCFrequency  (Hz)') 

legendOClassicalVSynthesis') 

change_model=inputCWould  you  like  to  change  your  model  again?  l=yes  2=no  '); 
end 
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fsynbase.m 


%  This  program  is  designed  to  calculate  FRFs  by  inverting  the  impedance  matrix. 

%  Changing  the  mass,  stiffness,  and  damping  matrices,  new  FRFs  will  be  calculated  by 

%  reformulating  the  matrices,  and  using  the  synthesis  method  in  the  frequency  domain. 

% 

clear 

%  loading  of  the  baseline  structure 

% 

disp('Select  the  file  which  contains  your  [M],  [K],  and  [C]  matrices' ) 

pause(2) 

[M_K_C,p]=uigetfile('*.mat,,,Load  [M],  [K],  [CD; 

load  (M_K_C) 

dofs=l:ndof; 

change_model=  1 ; 
while  change_model==l; 

changek=2;   changec=2;     changem=2;      %  initializes  logic  variables  to  'no'  status 

flag_k=0;    flag_c=0;     flag_m=0;    %  initializes  change  flag 

chk=0;    chc=0;     chm=0;     %  initializes  counter  for  #  of  changes  to  matrices 

cset=0;  %  initializes  cset  matrix  to  zero  to  avoid  null  index  error  in 

%  handling  the  spring  to  ground  if  no  changes  are  made  to 
%  the  structure. 

Hey_star=[];  HV_c=[];         %  (re)initializes  these  matrices  to  an  empty  matrix 

%  to  avoid  matrix  size  differences  each  time  the 
%  model  is  changed. 

K_c=K;    M_c=M;      C_c=C;  %  initializes  change  matrices  which  will  be 

%  used  in  the  classical  method,  to  the  original 
%  matrices. 

Kb=zeros(size(K,l));  %  initializes  the  spring,  and  damper  to  base  matrices 

Cb=zeros(size(C,l));  %  to  zero. 

%  Designation  where  base  excitation  is  located,  and  the  spring  constant  which 

%  connects  the  substructure  to  the  base  which  is  moving. 

% 

bset=input('At  what  dof(s)  are  your  structure  excited?  '); 

for  b_spring=l:length(bset); 

kb(b_spring)=input(fprintf ('input  the  value  of  the  spring  constant  which  connects  dof 

%g  to  the  base  *,bset(b_spring))); 
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end 

cb=zeros(  1  ,length(bset));  %  initializes  the  damper  constant  which  connects  the 

%  substructure  to  the  base  which  is  moving,  to  zero. 


%  Correcting  the  classical  [K]  &  [Kb]  matrices  to  include  the  spring  element 

%  connecting  the  substructure  to  the  base(s). 

% 

for  b_spring=l:length(bset); 

K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+... 

kb(b_spring); 
Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring))  +... 

kb(b_spring); 
end 

changek=inputCWould  you  like  to  change  your  stiffness  matrix?  l=yes  2=no  '); 
while  changek=l 

chk=chk+l ;  %  counter  for  determining  #  change  of  K  matrix 

flag_k=l; 

lk(chk)=input('input  value  of  added  spring  (lbs/in)  '); 

kdofl(chk)=input('input  constrained  dof  where  1st  end  of  added  spring  is  applied  '); 

kdof2(chk)=input('input  dof  where  2nd  end  of  added  spring  is  applied,  b  for  base,  or 

0  for  ground  '); 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 
%  are  not  included  in  the  bset.  A  change  from  dof  to  base  constitutes  an  addition  of 
%  excitation  points.  This  should  handled  when  asked  where  the  structure  is  excited 
% 

if  kdof  1  (chk)~=bset  &  kdof2(chk)=*b* 
while  kdofl(chk)~=bset  &  kdof2(chk)='b* 

disp('Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
dispCIf  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change=l 

error('Go  back  and  change  your  basic  model.') 
else 

kdof  l(chk)=input('input  constrained  dof  where  1st  end  of  added  damper  is 

applied  '); 
kdof2(chk)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b 

for  base,  or  0  for  ground  '); 
end 
end 
end 

%  updates  the  base  spring  constant  and  does  not  count  these  changes 
%  in  the  cset  matrix. 
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% 

if  (find(kdofl(chk)=bset))~=D  &  kdof2(chk)=='b' 
base=find(kdof  1  (chk)=bset); 
kb(base)  =  kb(base)  +  lk(chk); 
else 

if  kdof  l(chk)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  kdof  1  (chk)];  %  matrix  and  formation  of  new  c-set  matrix, 

end 
ifkdof2(chk)~=cset 

cset=[cset  kdof2(chk)]; 
end 
end 

%  Classical  method  of  reformulatting  the  [K]  matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk))  +  lk(chk); 

if  kdof2(chk)~=0  &  kdof2(chk)~='b' 

K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)4cdof2(chk))  -  lk(chk); 

K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)Jcdofl(chk))  -  lk(chk); 

K_c(kdof2(chk)Jcdof2(chk))=K_c(kdof2(chk)>kdof2(chk))  +  lk(chk); 
end 

%  Classical  method  of  reformulatting  [Kb]  matrix. 

% 

ifkdof2(chk)='b' 

kbb=lk(chk); 

Kb(kdofl(chk),kdofl(chk))=Kb(kdofl(chk),kdofl(chk))  +  kbb; 
end 

changek=input('Are  there  any  other  changes  to  the  stiffness  matrix?  l=yes  2=no  '); 
end 

changec=inputCWould  you  like  to  change  your  damping  matrix?  l=yes  2=no  '); 
while  changec=l 

chc=chc+l ;  %  counter  for  determining  #  changes  to  damping  matrix. 

flag_c=l; 

lc(chc)=input('input  value  of  added  damper  (lbs-secAn)  '); 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is  applied  '); 
cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b  for  base,  or 
0  for  ground  '); 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 
%  are  not  included  in  the  bset.  A  change  from  dof  to  base  constitutes  an  addition  of 
%  excitation  points.  This  should  handled  when  asked  where  the  structure  is  excited. 
% 

if  cdofl(chc)~=bset  &  cdof2(chc)==*b* 
while  cdofl(chc)~=bset  &  cdof2(chc)==*b' 

disp('Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
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model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change=l 

error('Go  back  and  change  your  basic  model.') 
else 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is 

applied  '); 
cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b 

for  base,  or  0  for  ground  '); 
end 
end 
end 


%  Updates  the  base  damper  constant  and  does  not  count  these  changes 

%  in  the  cset  matrix. 

% 

if  (find(cdofl(chc)==bset))~=[]  &  cdof2(chc)='b' 

base=find(cdof  1  (chc)=bset); 

cb(base)  =  cb(base)  +  lc(chc); 
else 

if  cdof  1  (chc)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  cdof  l(chc)];  %  matrix  and  formation  of  new  c-set  matrix. 

end 

if  cdof2(chc)~=cset 
cset=[cset  cdof2(chc)]; 

end 
end 

%  Classical  method  of  reformulatting  the  [C]  matrix. 

% 

C_c(cdof  1  (chc),cdof  1  (chc))=C_c(cdof  1  (chc),cdof  1  (chc))  +  lc(chc); 

if  cdof2(chc)~=0  &  cdof2(chc)~=*b' 

C_c(cdofl(chc),cdof2(chc))=C_c(cdofl(chc),cdof2(chc))  -  lc(chc); 

C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc))  -  lc(chc); 

C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc))  +  lc(chc); 
end 

%  Classical  method  of  reformulatting  [Cb]  matrix. 

% 

ifcdof2(chc)='b' 

cbb=lc(chc)* 

Cb(cdofl(chc),cdofl(chc))=Cb(cdofl(chc),cdofl(chc))  +  ebb; 
end 

changec=input('Are  there  any  other  changes  to  the  damping  matrix?  l=yes  2=no  '); 
end 

changem=inputCWould  you  like  to  change  your  mass  matrix?  l=yes  2=no  '); 
while  changem=l 
chm=chm+l; 
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flag_m=l; 

lm(chm)=input('input  value  of  added  mass  (lbs-secA2/in)  '); 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 

%  comparison  of  change  dof  with  the  c-set  matrix  and  formation  of 

%  new  c-set  matrix. 

% 

if  mdof (chm)~=cset 

cset=[cset  mdof(chm)]; 
end 

%  Classical  method  of  reformulatting  the  [M]  matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm))  +  lm(chm); 

changem=input('Are  there  any  other  changes  to  the  mass  matrix?  l=yes  2=no  '); 
end 

%  Handling  of  spring  added  to  ground 

% 

ground  =  find(cset=0); 

cset(ground)  =  Q; 

%  Formulation  of  iset  matrix  which  is  the  set  of  dofs  other  than 

%  those  in  the  cset  matrix  where  response  information  is  desired. 

% 

iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in?  '); 

%  Formation  of  zset  (c  U  b)  and  eset  (i  U  c  U  b)  matrices. 

% 

zset=[cset  bset];  eset=[iset  cset  bset]; 

%  Choosing  the  FRF  interested  in 

% 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

frf_index=find(eset==resp_dof) ; 

%  Choosing  the  type  of  input  and  out  FRF  interested  in 

% 

inp_type=input(ls  the  base  excitation  in  the  form  of:  l=displacement  or  2=acceleration 

'); 

resp_type=inputCWhat  type  of  response  are  you  interested  in:  l=displacement  or 

2=acceleration  '); 


%  Structural  Synthesis  method  of  calculating  FRF 


%  Calculation  of  unchanged  FRF 

% 

freq=.01:l:250+.01; 
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omega=2*pi*freq; 

j=length(omega); 

forw=l:j 

Z=K+i*omega(w)*C-omega(w)A2*M; 

H=inv(Z); 

%  Rearrangement  of  original  [H]  to  form  [Heb],  [Hez],  [Hzz],  [Hzb] 

% 

Heb=H(eset,bset); 

Hez=H(eset,zset) ; 

Hzz=H(zset,zset); 

Hzb=H(zset,bset); 

%  Formation  of  [del_Kc],  [del_Cc],  [del_Mc],  [del_Zc],  [del_Zb] 

%  &  [del_Z] 

% 

del_Kc=zeros(length(cset));         %  intializes  the  change  matrices 

del_Cc=zeros(length(cset));  %  to  zero  and  sets  their  sizes 

del_Mc=zerosQength(cset));         %  as  [cset  x  cset]  matrix 

ifflag_k=l 
fora=l:chk 

if  (find(kdofl(a)==bset))~=Q  &  kdof2(a)='b' 

indexkl=D;  %  provides  for  base  changes  to  not  be 

indexk2=Q;  %  included  in  del_Kc  matrix 

else 

indexk  1  =find(kdof  1  (a)=cset) ;  %  finds  where  in  cset  matrix 

indexk2=find(kdof2(a)==cset);  %  the  change  in  k  occurs 

end 

del_Kc(indexkl,  indexk  l)=del_Kc(indexkl,  indexk  1)  +  lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2)  -  lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl)  -  lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2)  +  lk(a); 
end 
end 

if  flag_c=l 
fora=l:chc 

if  (find(cdof l(a)==bset))~=D  &  cdof2(a)='b' 
indexcl=Q; 
indexc2=[]; 
else 

indexc  1  =find(cdof  1  (a)=cset);  %  finds  where  in  cset  matrix 

indexc2=find(cdof2(a)==cset);  %  the  change  in  C  occurs 

end 

del_Cc(indexcl, indexc  l)=del_Cc(indexcl, indexc  1)  +  lc(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexcl,indexc2)  -  lc(a); 
del_Cc(indexc2,indexcl)=del_Cc(indexc2,indexcl)  -  lc(a); 
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del_Cc(indexc2,indexc2)=del_Cc(indexc2,indexc2)  +  lc(a); 
end 
end 

if  flag_m=l 
fora=l:chm 
indexm=find(mdof(a)==cset);     %  finds  where  in  cset  matrix 

%  the  change  in  k  occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm)  +  lm(a); 
end 
end 

del_Zb=diag(kb+i*omega(w)*cb);  %  diag  is  used  in  the  event  >1  base 

%  excitation  is  present 

del_Zc=del_Kc+i*omega(w)*del_Cc-omega(w)A2*del_Mc; 

del_Z=[     del_Zc  zeros(size(del_Zc,l),size(del_Zb,2)) 

zeros(size(del_Zb,  1  ),size(del_Zc,2))         del_Zb] ; 


%  Base  excitation  Frequency  synthesis  equation 
% 

Hey_star(:,w)=(Heb*del_Zb-Hez*del_Z*inv(eye(size(Hzz,l))+Hzz*del_Z)*. 
Hzb*del_Zb)*ones(l  ,length(bset))*; 

if  inp_type==l  &  resp_type==2 

Hey_star(:,w)  =  Hey_star(:,w)*(-omega(w)A2); 
elseif  inp_type— 2  &  resp_type=l 

Hey_star(:,w)  =  Hey_star(:,w)*(-l/omega(w)A2); 
end 

%  Classical  method  of  calculating  the  FRF  by  reassembling  [Z]  and  taking  the 

%  inverse. 

% 

H=eye(size(M_c,  1 )); 

YV=ones(l,ndof)'; 

Z_c=(K_c+i*omega(w)*C_c-omega(w)A2*M_c); 

H_c=inv(Z_c); 

HV_c(:,w)  =  (H_c*(Kb+i*omega(w)*Cb))*YV; 

if  inp_type==l  &  resp_type=2 

HV_c(:,w)  =  HV_c(:,w)*(-omega(w)A2); 
elseif  inp_type=2  &  resp_type=l 

HV_c(:,w)  =  HV_c(:,w)*(-l/omega(w)A2); 
end 
end 
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%  Plotting  of  FRF  using  Classical  &  Frequency  Synthesis  Methods 

% 

semilogy(freq,abs(HV_c(resp_dof,:)),,r--',freq,abs(Hey_star(frf_index,:)),,b') 
titleflTrequency  Response  Function  (H',int2str(resp_dof),')'  ]) 
if  inp_type=l  &resp_type=2 
ylabelCFRF  Amplitude  (g/in)') 
elseif  inp_type=2  &  resp_type=l 

ylabel(*FRF  Amplitude  (in/g)*) 
else 

ylabelCFRF  Amplitude  (unitless)*) 
end 
xlabelCFrequency  (Hz)') 

legend('Classicar,'Syn  thesis') 

change_model=inputCWould  you  like  to  change  your  model  again?  l=yes  2=no  '); 
end 
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fsynmain.m 


%  This  program  is  designed  to  calculate  the  new  FRFs  of  a  system  after  changes  in  the 
%  mass,  stiffness,  and  damping  matrices  have  been  made.  The  new  FRFs  will  be 
%  calculated  using  the  synthesis  method  in  the  frequency  domain. 

% 

clear 

%  loads  the  baseline  structure  to  be  modified,  allows  for  modifications  to  the  baseline 

%  structure,  and  calls  the  function  delta.m 

% 

load  change 

%  Designation  of  the  frequency  range  and  step  size 

% 

%  initial      step  size      end 


freq_input  =  [     .01  1  250+.01  ]; 

freq=freq_input(l):freq_input(2):freq_input(3); 

omega=2*pi*freq; 

[wn,phi,zeta,Mmodal,phi_norm,Cmodal]  =  modal(M,K,C); 
[nee]  =  frfmodal(wn,phi,zeta,Mmodal,omega,eset); 

[h_star]  =  frfsynth(hee,kb,cb,omega,excite,eset,iset,bset,del_Kc,del_Cc,del_Mc); 

%  Choosing  the  type  of  input  and  output  FRF  interested  in 
% 

if  excite=2 
inp_type=input(ls  the  base  excitation  in  the  form  of:  l=displacement  or  2=acceleration 

'); 

resp_type=inputCWhat  type  of  response  are  you  interested  in:  l=displacement  or 

2=acceleration  '); 
if  inp_type=l  &resp_type=2 

h_star  =  h_star  *  (-omega(w)A2); 
elseif  inp_type=2  &  resp_type==l 

h_star  =  h_star  *  (-l/omega(w)A2); 
end 
end 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

if  excite=2 

inp_dof=[]; 
else 

inp_dof=inputCWhat  input  dof  are  you  interested  in?  '); 
end 
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[H_star_desired]  =  frf_sift(esetrresp_dof,inp_dof,excite,h_star); 


semilogy(freq',abs(H_star_desired),,g') 

tide(['Synthesized  Frequency  Response  Function  H',int2str([resp_dof  inp_dofj)  ]) 

ifinp_type=l  &resp_type=2 

ylabel(*FRF  Amplitude  (g/in)') 
elseif  inp_type==2  &  resp_type=l 

ylabelCFRF  Amplitude  (in/g)*) 
else 

ylabel(TFRF  Amplitude  (unidess)') 
end 
xlabel(Trequency  (rad/s)') 
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%  This  program  provides  for  the  loading  of  the  baseline  structure,  accepts  the  user  changes 
%  to  the  [M\,  [K],  and  [C]  matrices,  calls  the  function  delta.m,  and  stores  this  information 
%  to  be  loaded  from  another  program. 

% 

clear 


changek=2;   changec=2;     changem=2; 

status 

flag_k=0;    flag_c=0;     flag_m=0; 


%  initializes  logic  variables 


%  to  'no' 


chk=0;    chc=0;     chm=0; 


cset=0; 


bset=D; 


kdofl=D; 
cdofl=D; 
mdof=Q; 

lk=0;    lc=0; 

kbO=0; 


kdof2=D; 
cdof2=D; 


lm=0; 
cbO=0; 


b=num2str('b'); 


%  initializes  counter  for  #  of  changes  to  matrices 

%  initializes  cset  matrix  to  zero  to  avoid 
%  null  index  error  if  no  changes  are  made  to 
%  the  structure  and  bset  to  empty  matrix  to 
%  prevent  error  in  the  event  this  is  not  a 
%  base  excitation  and  bset  does  not  exist 

%  initalizes  the  matrices  which  keep 
%  track  of  the  dofs  where  changes  occur 


%  initializes  the  value  of  added  prarmeters 

%  initializes  kb  and  cb  to  zero  to  prevent  error  in  the  event 
%  this  is  not  a  base  excitation  and  kb  &  cb  does  not  exist 

%  allows  b  to  be  entered  as  a  numerical  value,  but 
%  declares  b  a  string  variable  to  used  for 
%  comparisons  in  logic  statements. 


disp('Are  there  already  reduced  FRF  &  IRF  matrices  in  the  correct  format  available,') 
disp('or  does  the  presynthesized  FRF  &IRF  matrices  need  to  be  generated  using  the  dof ') 
disp('locations  where  you  desire  to  change  the  structure?  ') 
FRF_IRF=input('l=reduced  FRF/IRF  already  exists        2=generate  reduced  FRF/IRF  '); 

%  Protects  against  user  making  an  error  in  choosing  how  FRF  &  IRF  is  obtained 
% 

ifFRF_IRF~=[12] 
while  FRF_IRF~=[  12] 
dispCError  in  choosing  how  FRF/IRF  is  obtained.  Choose  1  or  2.') 
FRF_IRF=input('l=reduced  FRF/IRF  already  exists  2=generate  reduced  FRF/IRF  '); 
end 
end 
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%  Loading  the  presynthesized  FRF/IRFs,  corresponing  frequency  vector,  and  vector  of 

%  dof  s  in  the  eset  if  reduced  FRF/IRF  already  exists.  Loading  of  M,  K,  &  C  matrices  if 

%  FRF/IRF  needs  to  be  generated 

% 

ifFRF_IRF=l 

dispCSelect  the  file  which  contains  your  presynthesized  FRF/IRFs,  a  corresponding 
frequency  ') 

disp(Vector,  and  the  vector  of  all  the  dofs  that  will  be  in  your  eset  vector.') 

pause(2) 
[hee_dofs_omega,p]=uigetfile('*.mat','Load  [FRF],  [IRF],  {dofs},  {omega}'); 

load  (hee_dofs_omega) 

else 

dispCSelect  the  file  which  contains  your  [M],  [K],  and  [C]  matrices' ) 

pause(2) 

[M_K_C,p]=uigetfile('*.mat','Load  [M],  [K],  [C]'); 

load  (M_K_C) 

dofs=l:ndof; 
end 

excite=inputCHow  is  the  structure  excited?  l=system    2=base  '); 

%  Protects  against  user  making  an  error  in  choosing  the  type  of  excitation. 
% 

if  excite~=[l  2] 
while  excite~=[l  2] 

dispCError  in  choosing  type  of  excitation.  Choose  1  or  2.') 
excite=inputCHow  is  the  structure  excited?  l=system    2=base  '); 
end 
end 

if  excite=2 
%  Designation  where  base  excitation  is  located. 
% 
bset=input('At  what  dof(s)  are  your  structure  excited?  '); 

%  Protects  against  user  making  an  error  in  choosing  the  location  of  the  base  excitation. 
% 

for  b_dof=l:length(bset) 
if  bset(b_dof)~=dof s 

while  bset(b_dof)~=dofs 

fprintfCError  in  choosing  location  of  base  excitation  at  dof  %g.\n',... 

bset(b_dof)) 
dispCEither  the  dof  chosen  does  not  exist  on  your  structure,  or  is  not ') 
disp('included  in  your  list  of  dofs  for  the  eset  vector.  ') 
bset(b_dof)=input(fprintf ('Choose  again  the  dof  for  base  exitation  %g  ',... 

b_dof)); 
end 
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end 
end 

%  Designation  of  the  spring  constant  which  connects  the  substructure  to  the  base  and 
%  protection  against  user  making  an  error  in  choosing  the  value  of  the  constant(s). 
% 

for  b_spring=l  :length(bset) 
kbO(b_spring)=input(fprintf(,input  the  value  of  the  spring  constant  which  connects  dof 

%g  to  the  base  \bset(b_spring))); 
kbO_zero=find(kbO=0); 
if  length(kbO)~=b_spring  I  kbO_zero~=Q 

while  length(kbO)~=b_spring  I  kbO_zero~=Q 

fprintf('You  made  an  error  in  entering  the  value  for  dof  %g  base 
excitation  spring  constant.\n',bset(b_spring)) 
kbO(b_spring)=0; 
kbO(b_spring)=input(fprintf(,Re-input  the  value  of  the  spring  constant 

which  connects  dof  %g  to  the  base  \bset(b_spring))); 
kbO_zero=find(kbO=0); 
end 
end 
end 

cbO=zeros(l  ,length(bset));  %  initializes  the  damper  constant  which  connects  the 

%  substructure  to  the  base  which  is  moving,  to  zero, 
end 

changek=inputeWould  you  like  to  make  stiffness  modifications  to  the  structure?  l=yes 

2=no  *); 
while  changek=l 

%  Designation  of  the  spring  change  and  protection  against  user  making  an  error  in 
%  choosing  the  value  of  the  constant(s). 
% 

chk=chk+l;     %  counter  for  determining  #  changes  to  stiffness  matrix. 
lk(chk)=input('input  value  of  added  spring  (lbs/in)  '); 
iflength(lk)~=chk 
while  length(lk)~=chk 

dispCYou  made  an  error  entering  the  value  for  the  stiffness  change.') 

lk(chk)=0; 

lk(chk)=input('Re-input  value  of  added  spring  (lbs/in)  '); 
end 
end 

%  User  selection  of  where  stiffness  changes  occur.  Protects  against  user  making  an 
%  error  in  choosing  the  locations  of  spring  changes. 
% 

kdofl(chk)=input('input  dof  where  1st  end  of  added  spring  is  applied  '); 
ifkdofl(chk)~=dofs 
while  kdofl(chk)~=dofs 

dispCError  in  choosing  location  of  1st  end  of  added  spring.  Either  the  dof 

chosen  does') 
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disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

kdofl(chk)=inputCRe-input  constrained  dof  where  1st  end  of  added  spring  is 

applied  '); 
end 
end 

kdof2(chk)=input('input  dof  where  2nd  end  of  added  spring  is  applied,  b  for  base,  or  0 

for  ground  '); 
if  kdof2(chk)~=dofs  &  kdof2(chk)~=0  &  kdof2(chk)~='b' 
while  kdof2(chk)~=dofs  &  kdof2(chk)~=0  &  kdof2(chk)~=*b' 

dispCError  in  choosing  location  of  2nd  end  of  added  spring.  Either  the  dof 

chosen  does') 
disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

kdof2(chk)=input('Re-input  dof  where  2nd  end  of  added  spring  is  applied,  b  for 

base,  or  0  for  ground  '); 
end 
end 

ifexcite=l  &  kdof2(chk)='b' 
while  excite==l  &  kdof2(chk)='b' 

dispCError  in  choosing  location  of  2nd  end  of  added  spring.  You  did  not 

designate  this') 
disp('as  a  base  excitation  structure.  Therefore  b  is  not  a  dof  option.') 
kdof2(chk)=input('Re-input  dof  where  2nd  end  of  added  spring  is  applied,  b  for 

base,  or  0  for  ground  '); 
end 
end 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 
%  are  not  included  in  the  bseL  A  change  from  dof  to  base  constitutes  an  addition  of 
%  excitation  points.  This  should  be  handled  when  asked  where  the  structure  is  excited. 
% 

if  excite==2  &  kdofl(chk)~=bset  &  kdof2(chk)='b' 
while  kdofl(chk)~=bset  &  kdof2(chk)='b' 

dispC'Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 
if  model_change=  1 

error('Go  back  and  change  your  basic  model.') 
else 
kdofl(chk)=input('Re-input  constrained  dof  where  1st  end  of  added  spring  is 

applied  '); 
kdof2(chk)=input('Re-input  dof  where  2nd  end  of  added  spring  is  applied,  b 

for  base,  or  0  for  ground  '); 
end 
end 
end 
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%  Provision  for  base  spring  constant  changes  to  not  be  included  in  the  cset  matrix. 

% 

if  excite=2  &  (find(kdofl(chk)==bset))~=[]  &  kdof2(chk)=='b' 

cset=cset; 
else 

if  kdof  1  (chk)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  kdof  1  (chk)] ;  %  matrix  and  formation  of  new  c-set  matrix. 

end 

ifkdof2(chk)~=cset 
cset=[cset  kdof2(chk)]; 

end 
end 

changek=input('Are  there  any  other  stiffness  modifications  to  the  structure?  l=yes  2=no'); 
end 

changec=inputCWould  you  like  to  make  damping  modifications  to  the  structure?  l=yes 

2=no  '); 

while  changec=l 

%  Designation  of  the  damper  change  and  protection  against  user  making  an  error  in 
%  choosing  the  value  of  the  constant(s). 
% 

chc=chc+l;  %  counter  for  determining  #  changes  to  damping  matrix. 

lc(chc)=input('input  value  of  added  damper  (lbs-sec/in)  '); 
if  length  (lc)~=chc 
while  length(lc)~=chc 

disp('You  made  an  error  entering  the  value  for  the  damper  change.') 

lc(chc)=0; 

lc(chc)=input('Re-input  value  of  added  damper  (lbs-sec/in)  '); 
end 
end 


%  User  selection  of  where  damper  changes  occur.  Protects  against  user  making  an  error 
%  in  choosing  the  locations  of  damper  changes. 
% 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is  applied  '); 
if  cdofl  (chc)~=dofs 
while  cdofl  (chc)~=dofs 

dispCError  in  choosing  location  of  1st  end  of  added  damper.  Either  the  dof 

chosen  does') 
disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

cdofl  (chc)=input('Re-input  constrained  dof  where  1st  end  of  added  damper  is 

applied  '); 
end 
end 

cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b  for  base,  or  0 

for  ground  '); 
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if  cdof2(chc)~=dofs  &  cdof2(chc)~=0  &  cdof2(chc)~=*b* 
while  cdof2(chc)~=dofs  &  cdof2(chc)~=0  &  cdof2(chc)~='b' 

dispCError  in  choosing  location  of  2nd  end  of  added  damper.  Either  the  dof 

chosen  does') 
disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dof  s  for  the') 
disp('eset  vector.  ' ) 

cdof2(chc)=inputCRe-input  dof  where  2nd  end  of  added  damper  is  applied,  b  for 

base,  or  0  for  ground  '); 
end 
end 

if  excite=l  &  cdof2(chc)=*b* 
while  excite==l  &  cdof2(chc)=='b' 

dispCError  in  choosing  location  of  2nd  end  of  added  damper.  You  did  not 

designate  this') 
disp('as  a  base  excitation  structure.  Therefore  b  is  not  a  dof  option.') 
cdof2(chc)=inputCRe-input  dof  where  2nd  end  of  added  damper  is  applied,  b  for 

base,  or  0  for  ground  '); 
end 
end 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 
%  are  not  included  in  the  bset.  A  change  from  these  dof  to  base  constitutes  an  addition 
%  of  excitation  points.  This  should  be  handled  when  asked  where  the  structure  is 
%  excited. 
% 

if  excite==2  &  cdofl(chc)~=bset  &  cdof2(chc)='b' 
while  cdofl(chc)~=bset  &  cdof2(chc)=='b* 

dispCConnecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change=  1 

error('Go  back  and  change  your  basic  model.') 
else 
cdof  l(chc)=input('Re-input  constrained  dof  where  1st  end  of  added  damper  is 

applied  '); 
cdof2(chc)=input('Re-input  dof  where  2nd  end  of  added  damper  is  applied,  b 

for  base,  or  0  for  ground  '); 
end 
end 
end 

%  Provision  for  base  damper  constant  changes  to  not  be  included  in  the  cset  matrix. 

% 

if  excite=2  &  (find(cdofl(chc)==bset))~=[]  &  cdof2(chc)='b, 

cset=cset; 
else 

if  cdof  1  (chc)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  cdof  l(chc)];  %  matrix  and  formation  of  new  c-set  matrix. 
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end 

if  cdof2(chc)~=cset 

cset=[cset  cdof2(chc)]; 
end 
end 

changec=input('Are  there  any  other  damping  modifications  to  the  structure?  l=yes  2=no'); 
end 

changem=inputCWould  you  like  to  make  a  mass  modification  to  the  structure?  l=yes 

2=no'); 

while  changem=l 

%  Designation  of  the  mass  change  and  protection  against  user  making  an  error  in 
choosing  die 

%  value  of  the  constant(s). 
% 

chm=chm+ 1 ;  %  counter  for  determining  #  changes  to  mass  matrix. 

lm(chm)=input('input  value  of  added  mass  (lbs-secA2/in)  '); 
if  length(lm)~=chm 
while  length(lm)~=chm 

disp('You  made  an  error  entering  the  value  for  the  mass  change.1) 
lm(chm)=0; 

lm(chm)=inputCRe-input  value  of  added  mass  (lbs-secA2/in)  '); 
end 
end 


%  User  selection  of  where  mass  changes  occur.  Protects  against  user  making  an  error  in 
%  choosing  the  locations  of  mass  changes. 
% 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 
if  mdof (chm)~=dof s 
while  mdof(chm)~=dofs 

dispCError  in  choosing  location  of  mass  change.  Either  the  dof  chosen  does 

not*) 
dispfexist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 
end 
end 

if  mdof(chm)~=cset  %  comparison  of  change  dof  with  the  c-set  matrix 

cset=[cset  mdof(chm)];  %  and  formation  of  new  c-set  matrix, 

end 

changem=input('Are  there  any  other  mass  modifications  to  the  structure?  l=yes  2=no  '); 
end 

%  Handling  of  spring  added  to  ground  to  eliminate  zero  from  cset 

% 
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ground  =  find(cset=0); 
cset(ground)  =  []; 

disp(The  following  is  a  list  of  dofs  were  you  have  made  changes  to  your  structure.') 

disp(This  represents  your  {cset}  vector:') 

disp(cset) 

%  Formulation  of  iset  matrix  which  is  the  set  of  dofs  other  than 
%  those  in  the  cset  matrix  where  response  information  is  desired. 
% 
ifFRF_IRF=l 

disp('Since  you  inputed  your  own  FRF/IRF  matrices,  your  iset  is  chosen  as  all  dofs 
remaining ') 

disp('in  your  eset/dof  vector  after  extracting  the  cset  vector.') 

iset=dofs; 

for  a=l  rlength(cset) 
extract(a)=find(cset(a)=dofs); 

end 

iset(extract)=Q; 
else 

iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in?  '); 
end 

%  Formation  of  zset  (c  U  b)  and  eset  (i  U  c  U  b)  matrices. 

% 

zset=[cset  bset];  eset=[iset  cset  bset]; 

[del_Kc,del_Cc,del_Mc,kb,cb]  = 

delta(cset,bset,kdof  1  ,kdof2,lk,cdof  1  ,cdof2,lc,mdof,lm,kb0,cb0); 

save  change  M  K  C  ndof  del_Kc  del_Cc  del_Mc  kb  cb  iset  cset  bset  zset  eset  excite 
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modal.m 


function  [wn,phi,zeta,Mmodal,phi_norm,Cmodal]  =  modal(M,K,C) 

%  This  function  is  designed  to  calculate  the  natural  frequencies  and 

%  mode  shapes 

% 

ndof=size(M,l); 

[PHI,lam]=eig(M\K); 

forj=l:ndof; 

lambda(j)=lam(j,j); 
end 

[wn,I]=sort(sqrt(lambda)); 
wn=rot90(wn,-l); 
forj=l:ndof; 

phi(:,j)=PHI(:,ia)); 
end 

%  Formulation  of  modal  mass  and  damping  matrices 

% 

Mmodal=phi'*M*phi; 

phi_norm^hi/(MmodalA.5); 

Mmodal=diag(Mmodal); 

Cmodal=phi'*C*phi; 

Cmodal=diag(Cmodal); 

%zetaO=0.0; 

%zeta=zetaO*ones(  1  ,length(Cmodal));  %  Constant  modal  damping  ratio 

zeta=Cmodal./(2*Mmodal.*wn);  %  Modal  damping  ratio  derived 

%  from  proportional  damping  matrix  [C] 
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frfmodal.m 


function  [h]  =  frfmodal(wn,phi,zeta,Mmodal,omega,firf_dof) 

%  This  function  is  designed  to  calculate  the  FRF  using  mode  shapes  and 
%  summing  the  FRF  over  a  number  of  modes. 

% 

ncol=aength(ftf_dof)*(length(frf_dof)+l))/2; 

%  Calculation  of  frequency  response  function 

% 

h=zeros(length(omega),ncol);  %  Initialization  to  zero,  and  sizing 


triangle  matrix. 


%  of  the  freq  x  lower 


h_old=0;  h_new=0;  %  initialization  of  the  matrices  to 

%  be  used  to  determine  when  modal 
%  convergence  has  occured. 

h_check=l;  %  Initialization  of  the  stop  criteria  for  summing  the 

%  modes 

nmodes=size(phi);       b=0;     %  Defining  the  number  of  modes  that  exist  and 

%  initializing  b  which  keeps  track  of  the 
%  of  modes  used  for  convergence  of  the  FRF. 

HOdprime_aa=zeros(length(frf_dof)); 

forb=l:nmodes 

%  Elimination  of  unnessecary  elements  and  rearrangement  of  original 

%  mode  shapes  {phi}  to  form  {phi_red}  and  [numjred]  which  is  now  an 

%  frf_dof  x  frf_dof  matrix. 

% 

num_^ed=phi(rrf_dof,b)*phi(f^f_dof,b),; 

for  w=l  ilength  (omega) 

den=(wn(b)A2  -  omega(w)A2  +  2*i*zeta(b)*wn(b)*omega(w))*Mmodal(b); 

H_red=num_red./den; 

%  Changes  the  [H_red]  FRFs  from  an  frf_dof  x  frf_dof  matrix  to 
%  a  freq  x  lower  triangle  matrix. 

% 

H_lowtri=tril(H_red);  %  pulls  out  the  lower  triangle 

%  and  places  zeros  everywhere  else 

symtry  =  find(H_lowtri==0);  %  finds  zero  values  in  the 
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%  lower  triangle  matrix 


H_lowtri(symtry)  =  []; 


h_mode(w,:)=H_lowtri; 


end 

h  =  h  +  h_mode; 
end 


%  deletes  all  values  of  zero  and 
%  turns  the  remaining  elements 
%  into  a  1  x  length(lower  triangle) 
%  vector 

%  saves  the  lower  triangle  vector 
%  at  each  freq.  These  FRFs  are  for 
%  a  single  mode  only. 


%  Summing  of  each  modal  FRF 
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frfsynth.m 


function  [h_star]  =  Msynth(heeJkb,cb,omega,excite,eset,iset,bset,... 
del_Kc,del_Cc,del_Mc) 

%  This  function  is  designed  to  calculate  the  new  FRF  using  the  synthesis 
%  method. 

% 

for  w=l  :length(omega) 
count=0;  %  initializes  the  counter  which  keeps  track  of 

%  which  column  of  the  hee  matrix  is  to  be 
%  placed  into  the  refomation  of  matrix  [Hee] 

for  col  =  lrlength(eset)  %  reforms  the  [Hee]  lower 

for  row  =  col:length(eset)  %  triangle  matrix 

count=count  +  1; 

Hee_lowtri(row,col)=hee(w,count); 
end 
end 

%  reforms  the  original  symmetric  [Hee]  matrix 

% 

Hee=Hee_lowtri; 

[sym_row,sym_col]  =  find(Hee_lowtri=0); 

for  n=l  :length(sym_row) 

Hee(sym_row(n),sym_col(n))=Hee(sym_col(n),sym_row(n)); 
end 

e=l:length(eset); 

c=length(iset)+ 1  :length(eset)-length(bset); 
b=length(eset)-length(bset)+ 1  :length(eset); 
z=length(iset)+l  :length(eset); 

%  Rearrangement  of  [Hee]  to  form  [Hee],  [Hcc],  [Hee]  if  the  structure 

%  is  system  excited. 

% 

if  excite=l 

Hec=Hee(e,c); 

Hce=Hee(c,e); 

Hcc=Hee(c,c); 
end 

%  Rearrangement  of  [Hee]  to  form  [Heb],  [Hez],  [Hzz],  [Hzb]  if  the 
%  structure  is  base  excited. 
% 

if  excite=2 
Heb=Hee(e,b); 
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Hez=Hee(e,z); 
Hzz=Hee(z,z); 
Hzb=Hee(z,b); 
end 


%  Formation  of  [del_Zc],  [del_Zb]  &  [del_Z] 

% 

del_Zc=del_Kc+i*omega(w)*del_Cc-omega(w)A2*del_Mc; 

if  excite=2 
del_Zb=diag(kb+i*omega(w)*cb);  %  diag  is  used  in  the  event 

%  >1  base  excitation 
%  is  present 

del_Z=[    del_Zc  zeros(size(del_Zc,l),size(del_Zb,2)) 

zeros(size(del_Zb,  1  ),size(del_Zc,2))         del_Zb] ; 
end 

%  Determining  which  synthesis  equations  are  to  be  used 

% 

if  excite==l 

%  Structure  excitation  Frequency  synthesis  equation 

% 

H_star=Hee-Hec*del_Zc*inv(eye(size(Hcc,l))  +  Hcc*del_Zc)*Hce; 

else 
%  Base  excitation  Frequency  synthesis  equation 
% 

h_star(:,w)=(Heb*del_Zb  - 
Hez*del_Z*inv(eye(size(Hzz,  1  ))+Hzz*del_Z)*Hzb*del_Zb)*ones(  1  ,length(bset))'; 

end 

if  excite=l 

%  Changes  the  [H*]  FRFs  from  an  eset  x  eset  matrix  to  a  freq  x  lower 
%  triangle  matrix,  and  creates  an  index  matrix  which  keeps  track  of 
%  the  original  coordinates  of  the  FRFs  in  the  [H*]  matrix. 

% 

Hstar_lowtri=tril(H_star);  %  pulls  out  the  lower  triangle  and  places 

%  zeros  everywhere  else 

symtry  =  find(Hstar_lowtri==0);  %  finds  zero  values  in  the  lower 

%  triangle  matrix 

Hstar_lowtri(symtry)  =  G;  %  deletes  all  values  of  zero  and  turns  the 
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%  remaining  elements  into  a  1  x  length(lower 
%  triangle)  vector 

h_star(w,:)=Hstar_lowtri;  %  saves  the  lower  triangle  vector  at  each  freq 


end 
end 


if  excite==l 

%  changes  h_star  matrix  so  freqs  are  in  each  column  and  the  FRFs  are 

%  in  the  rows 

% 

h_star=flipud(rot90(h_star)); 

end 
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frf  sift.m 


function  [H_star_desired]  =  rrf_sift(eset,resp_dof,inp_dof,excite,h_star) 


%  Locates  the  positions  in  the  eset  that  the  response  and  input  dofs  are  referring. 
%  These  locations  correspond  to  the  matrix  location  row  and  column  of  the  desired 
%  FRF  in  the  [H*]  matrix  which  is  eset  x  eset  and  partitioned.  Or  these  locations 
%  refer  to  the  row  of  the  [H*]  matrix  which  is  eset  x  1  for  a  base  excitation. 
%  These  indices  sift  through  the  [h*]  matrix  rows  to  find  the  desired  FRF. 
% 

%  Creates  an  index  matrix  which  keeps  track  of  the  original  coordinates  of  the  FRFs 

%  in  the  [H_star]  matrix. 

% 

count=0;  %  initializes  the  counter  which  keeps  track  of 

%  which  column  of  the  h  matrix  the  lower  triangle  [H_red] 

%  FRF  went  into 

for  col  =  1  :length(eset)  %  creates  a  matrix  which  matches  up 

for  row  =  cohlength(eset)  %  each  element  of  the  lower  triangle 

count=count  +  1;  %  H_star  matrix  with  the  rows  in  the 

hstar_index(count,:)=[row  col];  %  h_star  matrix  that  they  were  placed  in. 
end 
end 

frf_index=find(eset==resp_dof  I  eset==inp_dof); 

%  prevents  index  error  in  the  event  response  dof  and  input  dof  are  the  same 

% 

if  resp_dof=inp_dof 

frf_index=[frf_index  frfjndex]; 
end 

if  excite=l 

%  Uses  hstar_index  and  finds  the  row  in  the  lower  triangle  x  freq  matrix  that  the 
%  desired  FRF  is  located. 
% 

frf_row=find((frf_index(l)===hstar_index(:,l)  &  frf_index(2)=hstar_index(:,2))  I 
(frf_index(l)=hstar_index(:,2)  &  frf_index(2)=hstar_index(:,l))); 

H_star_desired  =  h_star(frf_row,:); 

else 
H_star_desired  =  h_star(frf_index,:); 

end 
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APPENDIX  B.  STATIC  DISPLACEMENT  SYNTHESIS  COMPUTER 

CODES 

The  following  is  a  list  and  a  brief  description  of  the  main  MATLAB 
computer  codes  that  were  written  in  order  to  perform  the  static  displacement 
synthesis  calculations. 

•  Guyan  xstat.m  -  uses  Guyan  reduction  methods  to  calculate  the  static 
displacements  on  a  structure. 

•  ssynmain.m  -  the  main  program  which  calls  the  modularized  static 
displacement  synthesis  programs  and  does  the  post  processing  of  the 
results. 

•  statchange.m  -  loads  the  baseline  structure  to  be  modified,  allows  for 
modifications  to  the  baseline  structure,  and  calls  the  function  statdelta.m 

•  statdelta.m  -  forms  the  modification  matrices  which  will  be  used  to  form 
impedance  matrices. 

•  modal. m  -  solves  for  the  structure's  natural  frequencies,  mode  shapes  and 
other  modal  matrices  and  vectors. 

•  frfmodal.m  -  calculates  the  structure's  FRFs. 

•  statsynth.m  -  performs  the  synthesis  on  the  baseline  structure  and  returns 

the  synthesized  static  displacements. 

The  full  codes  are  contained  on  subsequent  pages.  Any  codes  that  were 
previously  presented  will  not  be  repeated. 
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Guyanxstat.m 


%  The  purpose  in  this  program  is  to  use  Guyan  reduction  methods  in  order  to  calculate  the 

%  static  sisplacements  on  a  structure. 

% 

clear 

loadfe  add 


%  Stiffness  changes  to  the  structure 

% 

K(3,3)=K(3,3)-1000; 

K(93,93)=K(93,93)-1000; 

K(3,3)=K(3,3)+le4; 
K(93,93)=K(93,93)+le4; 

K(66,66)=K(66,66)+500; 

K(66,109)=K(66,109)-500; 

tic 

g=ones(size(M,l),l)*(-386.4); 

for  rot=  1:3:  size(M,  1  )-3 

g(rot)=0; 

end 

forrot=2:3:size(M,l)-2 

g(rot)=0; 
end 


K(18,18)=K(18,18)-1000; 
K(108,108)=K(108,108)-1000; 

K(18,18)=K(18,18)+le4; 
K(108,108)=K(108,108)+le4; 

K(109,109)=K(109,109)+500; 
K(109,66)=K(109,66)-500; 


F=M*g; 

dofset=l:ndof; 

oset=dofset; 

aset=[3  18  93  108  66  109]; 

for  a=l  :length(aset) 

kill(a)=find(aset(a)=dof set) ; 
end 

oset(kiU)=D; 
nset=[aset  oset]; 

Knn=K(nset,nset);      Fnn=F(nset); 
Kaa=K(aset,aset);       Kao=K(aset,oset); 
Toa=(- 1  )*inv(Koo)*Koa; 
T=[eye(length(aset));Toa]; 


Koo=K(oset,oset) ;      Koa=K(oset,aset) ; 


Kred_Guyan=T'*Knn*T;  Fred_Guyan=T'*Fnn; 

xstat_Guyan=Kred_Guyan\Fred_Guyan; 
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ssynmain.m 


%  This  program  is  designed  to  calculate  the  new  static  displacements  of  a  system  after 
%  changes  in  the  mass,  stiffness,  and  damping  matrices  have  been  made.  The  new  static 
%  displacements  will  be  calculated  using  the  static  displacement  synthesis  method 
% 
clear 

load  statchange 
stat_omega=.  1 :.  1  :.4; 
C=zeros(size(C)); 

[stat_wn,stat_phi,stat_zeta,stat_Mmodal]  =  modal(M,K,C); 

tO  =  clock; 

[stat_hee]  =  frfmodal(stat_wn,stat_phi,stat_zeta,stat_Mmodal,stat_omega,eset); 

[z_stat,Fred,Kred,Mred,HstarO_dp]  =  statsynth(stat_hee,kb,stat_omega,eset,iset,bset,.. 

del_Kc,del_Mc); 
tl  =  clock; 
stat_syn_time=etime(tl  ,t0) 
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statchange.m 


%  This  program  provides  for  the  loading  of  the  baseline  structure,  accepts  the  user  changes 
%  to  the  [M\,  [K],  and  [C]  matrices,  calls  the  function  statdelta.m,  and  stores  this 
%  information  to  be  loaded  from  another  program. 

% 

clear 


changek=2;   changec=2;     changem=2;     %  initializes  logic  variables 

status 

flag_k=0;    flag_c=0;     flag_m=0; 

chk=0;     chc=0;      chm=0; 


%  to  'no' 


cset=0; 


bset=D; 


kdofl=D; 
cdofl=[|; 
mdof=[]; 


kdof2=[|; 
cdof2=D; 


lk=0;    lc=0;     lm=0; 
kbO=0;  cb=0; 

b=num2str(*b*); 


%  initializes  counter  for  #  of  changes 
%  to  matrices 

%  initializes  cset  matrix  to  zero  to  avoid 
%  null  index  error  if  no  changes  are  made  to 
%  the  structure  and  bset  to  empty  matrix  to 
%  prevent  error  in  the  event  this  is  not  a 
%  base  excitation  and  bset  does  not  exist 

%  initalizes  the  matrices  which  keep 
%  track  of  the  dofs  where  changes  occur 


%  initializes  the  value  of  added  prarmeters 

%  initializes  kb  and  cb  to  zero  to 

%  prevent  error  in  the  event  this  is  not  a 

%  base  excitation  and  kb  &  cb  does  not  exist. 

%  allows  b  to  be  entered  as  a  numerical  value,  but 
%  declares  b  a  string  variable  to  used  for 
%  comparisons  in  logic  statements. 


%  Used  to  set  the  range  of  possible  dofs  for  the  user  to  choose  when  making  changes 

%  to  his  structure. 

% 

disp(ls  there  already  a  reduced  FRF  matrix  in  the  correct  format  available,1) 

disp('or  does  the  presynthesized  FRF  matrix  need  to  be  generated  using  the  dof ') 

disp('locations  where  you  desire  to  change  the  structure?  ') 

FRF=input('l=reduced  FRF  already  exists  2=generate  reduced  FRF  '); 

%  Protects  against  user  making  an  error  in  choosing  how  FRF  is  obtained. 
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% 

ifFRF~=[12] 
while  FRF~=[  12] 

dispCError  in  choosing  how  FRF  is  obtained.  Choose  1  or  2.') 
FRF=input('l=reduced  FRF  already  exists  2=generate  reduced  FRF  '); 

end 
end 

%  Loading  the  presynthesized  FRFs,  corresponing  frequency  vector,  and  vector  of  dofs  in 

%  the  eset  if  reduced  FRF  already  exists.  Loading  of  M,  K,  &  C  matrices  if  FRF  needs  to 

%  generated 

% 

ifFRF=l 

disp('Select  the  file  which  contains  your  presynthesized  FRFs,  a  corresponding 
frequency  ') 

disp('vector,  and  the  vector  of  all  the  dofs  that  will  be  in  your  eset  vector.') 

pause(2) 
[hee_dofs_omega,p]=uigetfile('*.mat7Load  [FRF],  {dofs},  {omega}'); 

load  (hee_dofs_omega) 

else 

disp('Select  the  file  which  contains  your  [M],  [K],  and  [C]  matrices' ) 

pause(2) 

[M_K_C,p]=uigetfile('*.mat','Load  [M],  [K],  [C]'); 

load  (M_K_C) 

dofs=l:ndof; 
end 

excite=inputCHow  is  the  structure  excited?  l=system    2=base  '); 

%  Protects  against  user  making  an  error  in  choosing  the  type  of  excitation. 
% 

if  excite~=[l  2] 
while  excite~=[l  2] 

dispCError  in  choosing  type  of  excitation.  Choose  1  or  2.') 
excite=inputCHow  is  the  structure  excited?  l=system    2=base  '); 
end 
end 

if  excite=2 

%  Designation  where  base  excitation  is  located 

% 

bset=input('At  what  dof(s)  are  your  structure  excited?  '); 
%    bset=l; 

%  Protects  against  user  making  an  error  in  choosing  the  location  of  the  base  excitation. 
% 

for  b_dof=l  :length(bset) 
if  bset(b_dof)~=dofs 

while  bset(b_dof)~=dofs 

fprintfCError  in  choosing  location  of  base  excitation  at  dof  %g.\n',... 
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bset(b_dof)) 
dispCEither  the  dof  chosen  does  not  exist  on  your  structure,  or  is  not ') 
disp('included  in  your  list  of  dofs  for  the  eset  vector.  ') 
bsetCb.dof^nputCfprintfOChoose  again  the  dof  for  base  exitation  %g  ',... 

b_dof)); 
end 
end 
end 

%  Designation  of  the  spring  constant  which  connects  the  substructure  to  the  base  which 
is  moving. 
%  and  protection  against  user  making  an  error  in  choosing  the  value  of  the  constant(s). 
% 

for  b_spring=l:length(bset) 
kbO(b_spring)=input(fprintf('input  the  value  of  the  spring  constant  which  connects  dof 

%g  to  the  base  \bset(b_spring))); 
kbO_zero=find(kbO=0); 
if  length(kbO)~=b_spring  I  kbO_zero~=n 

while  length(kbO)~=b_spring  I  kbO_zero~=[] 

fprintf('You  made  an  error  in  entering  the  value  for  dof  %g  base 
excitation  spring  constant.\n',bset(b_spring)) 
kbO(b_spring)=0; 
kbO(b_spring)=input(fprintf('Re-input  the  value  of  the  spring  constant 

which  connects  dof  %g  to  the  base  \bset(b_spring))); 
kbO_zero=find(kbO=0); 
end 
end 
end 

cb=zeros(l  ,length(bset));  %  initializes  the  damper  constant  which  connects  the 

%  substructure  to  the  base  which  is  moving,  to  zero, 
end 

changek=inputCWould  you  like  to  make  stiffness  modifications  to  the  structure?  l=yes 

2=no  '); 

while  changek=l 

%  Designation  of  the  spring  change  and  protection  against  user  making  an  error  in 
choosing  the 

%  value  of  the  constant(s). 
% 

chk=chk+l;     %  counter  for  determining  #  changes  to  stiffness  matrix. 
lk(chk)=input('input  value  of  added  spring  (lbs/in)  '); 
if  length(lk)~=chk 
while  length(lk)~=chk 

dispCYou  made  an  error  entering  the  value  for  the  stiffness  change.') 
lk(chk)=0; 

lk(chk)=input('Re-input  value  of  added  spring  (lbs/in)  '); 
end 
end 
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%  User  selection  of  where  stiffness  changes  occur.  Protects  against  user  making  an 
%error  in  choosing 
%  the  locations  of  spring  changes. 
% 

kdof  l(chk)=input('input  constrained  dof  where  1st  end  of  added  spring  is  applied  '); 
ifkdofl(chk)~=dofs 
while  kdofl(chk)~=dofs 

dispCError  in  choosing  location  of  1st  end  of  added  spring.  Either  the  dof 

chosen  does') 
disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

kdofl(chk)=inputCRe-input  constrained  dof  where  1st  end  of  added  spring  is 

applied  '); 
end 
end 

kdof2(chk)=input('input  dof  where  2nd  end  of  added  spring  is  applied,  b  for  base,  or  0 
for  ground  '); 

if  kdof2(chk)~=dofs  &  kdof2(chk)~=0  &  kdof2(chk)~='b' 
while  kdof2(chk)~=dofs  &  kdof2(chk)~=0  &  kdof2(chk)~=*b* 

dispCError  in  choosing  location  of  2nd  end  of  added  spring.  Either  the  dof 

chosen  does') 
disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

kdof2(chk)=input('Re-input  dof  where  2nd  end  of  added  spring  is  applied,  b  for 

base,  or  0  for  ground  '); 
end 
end 

ifexcite=l  &  kdof2(chk)==*b' 
while  excite==l  &  kdof2(chk)=='b' 

dispCError  in  choosing  location  of  2nd  end  of  added  spring.  You  did  not 

designate  this') 
disp('as  a  base  excitation  structure.  Therefore  b  is  not  a  dof  option.') 
kdof2(chk)=input('Re-input  dof  where  2nd  end  of  added  spring  is  applied,  b  for 

base,  or  0  for  ground  '); 
end 
end 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 
%  are  not  included  in  the  bseL  A  change  from  dof  to  base  constitutes  an  addition  of 
%  excitation  points.  This  should  be  handled  when  asked  where  the  structure  is  excited. 
% 

if  excite==2  &  kdof  1  (chk)~=bset  &  kdof2(chk)=='b' 
while  kdofl(chk)~=bset  &  kdof2(chk)=='b' 

dispCConnecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=input('Was  this  change  correct?  l=yes  2=no  '); 

if  model_change=l 
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error('Go  back  and  change  your  basic  model.') 
else 
kdofl(chk)=input('Re-input  constrained  dof  where  1st  end  of  added  spring  is 

applied  '); 
kdof2(chk)=input('Re-input  dof  where  2nd  end  of  added  spring  is  applied,  b 

for  base,  or  0  for  ground  '); 
end 
end 
end 

%  Provision  for  base  spring  constant  changes  to  not  be  included  in  the  cset  matrix. 

% 

if  excite=2  &  (find(kdofl(chk)==bset))~=[]  &  kdof2(chk)=='b' 

cset=cset; 
else 

if  kdof  1  (chk)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  kdof  l(chk)];  %  matrix  and  formation  of  new  c-set  matrix. 

end 

ifkdof2(chk)~=cset 
cset=[cset  kdof2(chk)]; 

end 
end 

changek=input('Are  there  any  other  stiffness  modifications  to  the  structure?  l=yes  2=no 

'); 

end 

changec=dnputCWould  you  like  to  make  damping  modifications  to  the  structure?  l=yes 

2=no  *); 

while  changec=l 

%  Designation  of  the  damper  change  and  protection  against  user  making  an  error  in 
choosing  the 

%  value  of  the  constant(s). 
% 

chc=chc+l;  %  counter  for  determining  #  changes  to  damping  matrix. 

lc(chc)=input('input  value  of  added  damper  (lbs-sec/in)  '); 
if  length(lc)~=chc 
while  length(lc)~=chc 

disp('You  made  an  error  entering  the  value  for  the  damper  change.') 
lc(chc)=0; 

lc(chc)=input('Re-input  value  of  added  damper  (lbs-sec/in)  '); 
end 
end 


%  User  selection  of  where  damper  changes  occur.  Protects  against  user  making  an  error 

%  in  choosing  the  locations  of  damper  changes. 

% 

cdofl(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is  applied  '); 

if  cdof  1  (chc)~=dofs 
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while  cdofl(chc)~=dofs 

dispCError  in  choosing  location  of  1st  end  of  added  damper.  Either  the  dof 
chosen  does') 

disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dof  s  for  the') 
disp('eset  vector.  ' ) 

cdof  l(chc)=input('Re-input  constrained  dof  where  1st  end  of  added  damper  is 
applied  '); 
end 
end 

cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b  for  base,  or  0 
for  ground  '); 

if  cdof2(chc)~=dofs  &  cdof2(chc)~=0  &  cdof2(chc)~=*b' 
while  cdof2(chc)~=dofs  &  cdof2(chc)~=0  &  cdof2(chc)~='b' 

dispCError  in  choosing  location  of  2nd  end  of  added  damper.  Either  the  dof 
chosen  does') 

disp('not  exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
disp('eset  vector.  ' ) 

cdof2(chc)=input('Re-input  dof  where  2nd  end  of  added  damper  is  applied,  b  for 
base,  or  0  for  ground  '); 
end 
end 

if  excite=l  &  cdof2(chc)='b' 
while  excite==l  &  cdof2(chc)=='b' 

dispCError  in  choosing  location  of  2nd  end  of  added  damper.  You  did  not 
designate  this') 

disp('as  a  base  excitation  structure.  Therefore  b  is  not  a  dof  option.') 
cdof2(chc)=input('Re-input  dof  where  2nd  end  of  added  damper  is  applied,  b  for 
base,  or  0  for  ground  '); 
end 
end 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 
%  are  not  included  in  the  bset.  A  change  from  these  dof  to  base  constitutes  an  addition 
of 

%  excitation  points.  This  should  be  handled  when  asked  where  the  structure  is  excited. 
% 

if  excite=2  &  cdofl(chc)~=bset  &  cdof2(chc)='b' 
while  cdofl(chc)~=bset  &  cdof2(chc)=='b* 

dispCConnecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 
point.') 

model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change==l 

error('Go  back  and  change  your  basic  model.') 
else 
cdof  l(chc)=input('Re-input  constrained  dof  where  1st  end  of  added  damper  is 
applied  '); 
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cdof2(chc)=input('Re-input  dof  where  2nd  end  of  added  damper  is  applied,  b 
for  base,  or  0  for  ground  '); 
end 
end 
end 

%  Provision  for  base  damper  constant  changes  to  not  be  included  in  the  cset  matrix. 

% 

if  excite=2  &  (find(cdofl(chc)==bset))~=G  &  cdof2(chc)=V 

cset=cset; 
else 

if  cdof  1  (chc)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  cdof  l(chc)];  %  matrix  and  formation  of  new  c-set  matrix. 

end 

if  cdof2(chc)~=cset 

cset=[cset  cdof2(chc)]; 

end 
end 

changec=input('Are  there  any  other  damping  modifications  to  the  structure?  l=yes  2=no 

end 

changem=inputCWould  you  like  to  make  a  mass  modification  to  the  structure?  l=yes  2=no 

while  changem=l 

%  Designation  of  the  mass  change  and  protection  against  user  making  an  error  in 
choosing  the 

%  value  of  the  constant(s). 
% 

chm=chm+ 1 ;  %  counter  for  determining  #  changes  to  mass  matrix. 

lm(chm)=input('input  value  of  added  mass  (lbs-secA2/in)  '); 
if  length(lm)~=chm 
while  length(lm)~=chm 

disp('You  made  an  error  entering  the  value  for  the  mass  change.') 
lm(chm)=0; 

lm(chm)=input('Re-input  value  of  added  mass  (lbs-secA2/in)  '); 
end 
end 


%  User  selection  of  where  mass  changes  occur.  Protects  against  user  making  an  error  in 
%  choosing  the  locations  of  mass  changes. 
% 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 
if  mdof(chm)~=dof s 
while  mdof  (chm)~=dofs 

dispCError  in  choosing  location  of  mass  change.  Either  the  dof  chosen  does 
not') 

disp('exist  on  your  structure,  or  is  not  included  in  your  list  of  dofs  for  the') 
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disp('eset  vector.  ' ) 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 
end 
end 

if  mdof(chm)~=cset  %  comparison  of  change  dof  with  the  c-set  matrix 

cset=[cset  mdof(chm)];  %  and  formation  of  new  c-set  matrix, 

end 

changem=input('Are  there  any  other  mass  modifications  to  the  structure?  l=yes  2=no  '); 
end 

%  Handling  of  spring  added  to  ground  to  eliminate  zero  from  cset 

% 

ground  =  find(cset=0); 

cset(ground)  =  []; 

disp(The  following  is  a  list  of  dofs  were  you  have  made  changes  to  your  structure.') 

disp(This  represents  your  {cset}  vector:') 

disp(cset) 

%  Formulation  of  iset  matrix  which  is  the  set  of  dofs  other  than 
%  those  in  the  cset  matrix  where  response  information  is  desired. 
% 

ifFRF=l 
disp('Since  you  inputed  your  own  FRF  matrix,  your  iset  is  chosen  as  all  dofs  remaining 

") 
disp('in  your  eset/dof  vector  after  extracting  the  cset  vector.') 

iset=dofs; 

for  a=l:length(cset) 

extract(a)=find(cset(a)==dofs); 

end 

iset(extract)=n; 
else 

iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in?  '); 
end 

%  Formation  of  zset  (c  U  b)  and  eset  (i  U  c  U  b)  matrices. 

% 

zset=[cset  bset];  eset=[iset  cset  bset]; 

[del_Kc,dd_Mc,kb]  =  statdelta(cset,bset,kdofl,kdof2,lk,mdof,lm,kb0); 
save  statchange  MKC  ndof  del_Kc  del_Mc  kb  iset  cset  bset  zset  eset  excite 
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statdelta.m 


%  The  purpose  of  this  function  is  to  form  the  modification  matrices 

% 

function  [del_Kc,del_Mc,kb]  =  statdelta(cset,bset,kdofljcdof2,lk,mdof,lmjcb0) 

%  Formation  of  [del_Kc],  [del_Cc],  &  [del_Mc]  and  updating  of  kb  &  cb 

% 


del_Kc=zeros(length(cset)); 
del_Mc=zeros(length(cset)); 

matrix 

kb=kbO; 


%  intializes  the  change  matrices 
%  to  zero  and  sets  their  si2es 

%  as  [cset  x  cset] 


%  initializes  the  base  excitation  spring  constant 
%  to  the  initial  value  inputed  when  structure  was 
%  formed. 


ifkdofl~=D 
for  a=l:length(kdofl) 

if  (find(kdofl(a)=bset))~=Q  &  kdof2(a)=='b' 
base=find(kdof  1  (a)=bset); 
kb(base)  =  kb(base)  +  lk(a); 


else 


end 
end 
end 


indexkl=find(kdofl(a)==cset);    %  finds  where  in  cset  matrix 
indexk2=find(kdof2(a)=cset);    %  the  change  in  k  occurs 

del_Kc(indexkl,indexkl)=del_Kc(indexkl,indexkl)  +  lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2)  -  lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl)  -  lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2)  +  lk(a); 


ifmdof~=[] 
for  a=l:length(mdof) 

indexm=find(mdof (a)==cset) ; 


%  finds  where  in  cset  matrix 
%  the  change  in  M  occurs 


del_Mc(indexm,indexm)=del_Mc(indexm,indexm)  +  lm(a); 
end 
end 
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statsynth.m 


function  [z_stat,Fred,Kred,Mred,HstarO_dp]  =  statsynth(hee,kb,omega,eset,iset,bset,. 

del_Kc,del_Mc) 

%  This  function  is  designed  to  calculate  the  new  static  displacement  using  the  synthesis 
%  method. 

% 

puii=D; 

for  w=l:length(omega) 
count=0;  %  initializes  the  counter  which  keeps  track  of 

%  which  column  of  the  hee  matrix  is  to  be 
%  placed  into  the  refomation  of  matrix  [Hee] 

for  col  =  l:length(eset)  %  reforms  the  [Hee]  lower 

for  row  =  col:length(eset)  %  triangle  matrix 

count=count  +  1; 

Hee_lowtri(row,col)=hee(w,count); 
end 
end 

%  reforms  the  original  symmetric  [Hee]  matrix 

% 

Hee=Hee_lowtri; 

[sym_row,sym_col]  =  find(Hee_lowtri=0); 

for  n=l  :length(sym_row) 

Hee(sym_row(n),sym_col(n))=Hee(sym_col(n),sym_row(n)); 
end 

e=l:length(eset); 

c=length(iset)+ 1  :length(eset);  %-length(bset); 
%    b=length(eset)-length(bset)+ 1  rlength(eset); 
%    z=length(iset)+ 1  :length(eset); 

%  Rearrangement  of  [Hee]  to  form  [Hee],  [Hcc],  [Hee]  if  the  structure 
%  is  system  excited. 
% 
%    ifexcite=l 
Hec=Hee(e,c); 
Hce=Hee(c,e); 
Hcc=Hee(c,c); 
%    end 

%  Formation  of  [del_Zc],  [del_Zb]  &  [del_Z] 

% 
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del_Zc=del_Kc-omega(w)A2*del_Mc; 

del_Zb=diag(kb);  %  diag  is  used  in  the  event      >1  base  excitation 

%  is  present 

del_Z=[    del_Zc  zeros(size(del_Zc,l),size(del_Zb,2)) 

zeros(size(del_Zb,  1 )  ,size(del_Zc  ,2))         delJZb] ; 


%  Structure  excitation  Frequency  synthesis  equation 

% 

H_star=Hee-Hec*del_Z*inv(eye(size(Hcc,l))  +  Hcc*del_Z)*Hce; 

%  Checks  to  see  if  there  are  any  redundant  dofs  in  the  eset  due  to  changes  and 

%  bset  occurring  at  the  same  dofs.  If  there  is  redundancy,  it  is  located  in 

%  eset,  and  these  rows  and  columns  are  deleted  from  H_star.  This  is  done 

%  because  the  redundant  dof  column(s)  and  row(s)  in  H_star  creates  a  singular 

%  matrix  which  is  not  invertable  and  Kred  cannot  be  determined 

% 

for  u=l:length(bset) 

locate=find(bset(u)=eset); 

if  length(locate)  >  1 
pull=[pull  locate(2)]; 

end 
end 

H_star(pull,:)=n; 
H_star(:,pull)=D; 

ifw=l 

Kred=inv(H_star); 

H0=H_star; 
end 
ifw=2 

Hl=H_star, 
end 
ifw=3 

H2=H_star; 
end 
if  w=4 

H3=H_star; 
end 

Hstar0_dp=(-H3+4*H2-5*Hl+2*H0)./(omega(2)-omega(l))A2; 
Mred=.5*Kred*HstarO_dp*Kred; 
end 

ga=ones(size(Mred,l ),  1  )*(-386.4); 

Fred=Mred*ga; 

z_stat=KrecNFred; 
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APPENDIX  C  TIME  DOMAIN  SYNTHESIS  COMPUTER  CODES 

The  following  is  a  list  and  a  brief  description  of  the  main  MATLAB 
computer  codes  that  were  written  in  order  to  perform  the  time  domain 
synthesis  calculations. 

•  tsynconv.m  -  performs  the  time  domain  synthesis  on  a  structure  which 
experiences  a  base  excitation.  The  dynamic  responses  are  compared  to  the 
responses  calculated  using  the  classical  convolution  integral  method. 

•  tsynode45.m  -  performs  the  time  domain  synthesis  on  a  structure  which 
experiences  a  base  excitation.  The  dynamic  responses  are  compared  to  the 
responses  calculated  using  the  classical  direct  integration  method. 

•  tsynmain.m  -  the  main  program  which  calls  the  modularized  time 
synthesis  programs  and  does  the  post  processing  of  the  results. 

•  change. m  -  loads  the  baseline  structure  to  be  modified,  allows  for 
modifications  to  the  baseline  structure,  and  calls  the  function  delta.m 

•  delta.m  -  forms  the  modification  matrices  which  will  be  used  to  form 
impedance  matrices  or  determine  coupling  forces. 

•  modal.m  -  solves  for  the  structure's  natural  frequencies,  mode  shapes  and 
other  modal  matrices  and  vectors. 

•  impmodal.m  -  calculates  the  structure's  IRFs. 

•  buildA.m  -  builds  the  quadrature  matrices. 
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•  fBlastForcing.m  -  inputs  the  base  excitation  as  a  function  of  time. 

•  timesynth.m  -  performs  the  synthesis  on  the  baseline  structure  and 
returns  the  synthesized  transient  responses. 

•  time  sift.m  -  selects  the  dynamic  response  that  the  user  wants  to  perform 

post  processing  on. 

The  full  codes  are  contained  on  subsequent  pages.  Any  codes  that  were 
previously  presented  will  not  be  repeated. 
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tsynconv.m 


%  The  purpose  of  this  program  is  to  perform  the  time  domain  synthesis  on  a  structure  and 
%  compare  the  dynamic  responses  to  those  calculated  using  a  convolution  integral  method 

% 

clear 

disp('Select  the  file  which  contains  your  [M],  [K],  and  [C]  matrices' ) 
pause(2) 

[M_K_C,p]=uigetfile(,*.mat',*Load  [M],  [K],  [C]'); 
load  (M_K_C) 
dofs=l:ndof; 
C=00*C; 

[wn,phi,zeta,Mmodal,phi_norm]  =  modal(M,K,C); 

changek=2;   changec=2;     changem=2;     %  initializes  logic  variables 

%  to  'no'  status 


flag_k=0;    flag_c=0;     flag_m=0; 
chk=0;     chc=0;      chm=0; 


%  initializes  counter  for  #  of  changes 
%  to  matrices 


cset=0; 


%  initializes  cset  matrix  to  zero  to  avoid  null  index  error  in 
%  handling  the  spring  to  ground  if  no  changes  are  made  to 
%  the  structure. 


Hey_star=D; 


HV_c=[]; 


K_c=K;     M_c=M;      C_c=C; 

Kb=zeros(size(K,  1 )); 
Cb=zeros(size(C,l)); 

K_vec  =  zeros(ndof,l); 

b=num2str(,b'); 

clearMKC 


%  (re)initializes  these  matrices  to  an  empty  matrix 
%  to  avoid  matrix  size  differences  each  time  the 
%  model  is  changed. 

%  initializes  change  matrices  which  will  be  used  in 
%  the  classical  method,  to  the  original  matrices. 

%  initializes  the  spring,  and  damper  to  base  matrices 
%  to  zero. 


%  allows  b  to  be  entered  as  a  numerical  value,  but 
%  declares  b  a  string  variable  to  used  for 
%  comparisons  in  logic  statements. 


%  Designation  where  base  excitation  is  located,  and  the  spring  constant  which 
%  connects  the  substructure  to  the  base  which  is  moving. 
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% 

bset=input('At  what  dof(s)  are  your  structure  excited?  '); 
for  b_spring=l:length(bset); 
kb(b_spring)=input(fprintf('input  the  value  of  the  spring  constant  which  connects  dof 

%g  to  the  base  ',bset(b_spring))); 
end 
cb=zeros(l  ,length(bset));  %  initializes  the  damper  constant  which  connects  the 

%  substructure  to  the  base  which  is  moving,  to  zero. 

%  Correcting  the  classical  [K]  &  [Kb]  matrices  to  include  the  spring  element 

%  connecting  the  substructure  to  the  base(s). 

% 

for  b_spring=l:length(bset); 

K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+kb(b_spring); 

Kb(bset(D_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring))  + 

kb(b_spring); 

K_vec(bset(b_spring))=K_vec(bset(b_spring))  +  kb(b_spring); 
end 

changek=inputCWould  you  like  to  change  your  stiffness  matrix?  l=yes  2=no  '); 
while  changek=l 

chk=chk+ 1 ;  %  counter  for  determining  #  changes 

flag_k=l; 

lk(chk)=input('input  value  of  added  spring  (lbs/in)  '); 

kdofl(chk)=input('input  constrained  dof  where  1st  end  of  added  spring  is  applied  '); 

kdof2(chk)=input('input  dof  where  2nd  end  of  added  spring  is  applied,  b  for  base,  or 

0  for  ground  '); 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 

%  are  not  included  in  the  bset.  A  change  from  dof  to  base  constitutes  an  addition  of 

%  excitation  points.  This  should  handled  when  asked  where  the  structure  is  excited. 

% 

if  kdofl(chk)~=bset  &  kdof2(chk)='b' 

while  kdofl(chk)~=bset  &  kdof2(chk)=='b' 
disp('Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
dispCIf  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change=l 

error('Go  back  and  change  your  basic  model.') 
else 

kdofl(chk)=input('input  constrained  dof  where  1st  end  of  added 

damper  is  applied  '); 
kdof2(chk)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b 
for  base,  or  0  for  ground  '); 
end 
end 
end 
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%  updates  the  base  spring  constant  and  does  not  count  these  changes 

%  in  the  cset  matrix. 

% 

if  (find(kdofl(chk)=bset))~=Q  &  kdof2(chk)=='b' 

base=find(kdof  1  (chk)=bset); 

kb(base)  =  kb(base)  +  lk(chk); 

K_vec(kdofl(chk))=K_vec(kdofl(chk))  +  lk(chk); 
else 
if  kdof  1  (chk)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  kdofl(chk)];  %  matrix  and  formation  of  new  c-set  matrix. 

end 

if  kdof2(chk)~=cset 
cset=[cset  kdof2(chk)]; 

end 
end 

%  Classical  method  of  reformulatting  the  [K]  matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk))  +  lk(chk); 

if  kdof2(chk)~=0  &  kdof2(chk)~='b' 

K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)Jcdof2(chk))  -  lk(chk); 

K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)Jcdofl(chk))  -  lk(chk); 

K_c(kdof2(chk)4cdof2(chk))=K_c(kdof2(chk)^dof2(chk))  +  lk(chk); 
end 

%  Classical  method  of  reformulatting  [Kb]  matrix. 

% 

ifkdoOCchk^'b* 

kbb=lk(chk); 

Kb(kdofl(chk),kdof  l(chk))=Kb(kdofl(chk),kdofl(chk))  +  kbb; 
end 

changek=input('Are  there  any  other  changes  to  the  stiffness  matrix?  l=yes  2=no  '); 

end 

changec=inputCWould  you  like  to  change  your  damping  matrix?  l=yes  2=no  '); 
while  changec=l 

chc=chc+ 1 ;  %  counter  for  determining  #  columns  in  mapping  matrix. 

flag_c=l;        %  mapping  matrix  flag 

lc(chc)=input('input  value  of  added  damper  (lbs-sec/in)  '); 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is  applied  '); 
cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b  for  base,  or 
0  for  ground  '); 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 

%  are  not  included  in  the  bset.  A  change  from  dof  to  base  constitutes  an  addition  of 

%  excitation  points.  This  should  handled  when  asked  where  the  structure  is  excited. 

% 

if  cdofl(chc)~=bset  &  cdof2(chc)=='b' 
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while  cdofl(chc)~=bset  &  cdof2(chc)='b' 
disp('Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 
point.') 

model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change=l 

error('Go  back  and  change  your  basic  model.') 
else 

cdofl(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is 

applied  '); 
cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b 

for  base,  or  0  for  ground  '); 
end 
end 
end 

%  Updates  the  base  damper  constant  and  does  not  count  these  changes 

%  in  the  cset  matrix. 

% 

if  (find(cdofl(chc)=bset))~=[]  &  cdof2(chc)='b' 

base=find(cdof  1  (chc)=bset); 

cb(base)  =  cb(base)  +  lc(chc); 

C_vec(cdofl(chc))=C_vec(cdofl(chc))  +  lc(chc); 


else 


if  cdof  1  (chc)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  cdof  1  (chc)] ;        %  matrix  and  formation  of  new  c-set  matrix, 
end 

if  cdof2(chc)~=cset 
cset=[cset  cdof2(chc)]; 
end 


end 


%  Classical  method  of  reformulatting  the  [C]  matrix. 

% 

C_c(cdof  1  (chc),cdof  1  (chc))=C_c(cdof  1  (chc),cdof  1  (chc))  +  lc(chc); 

if  cdof2(chc)~=0  &  cdof2(chc)~='b* 

C_c(cdofl(chc),cdof2(chc))=C_c(cdofl(chc),cdof2(chc))  -  lc(chc); 

C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc))  -  lc(chc); 

C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc))  +  lc(chc); 
end 

%  Classical  method  of  reformulating  [Cb]  matrix. 

% 

ifcdof2(chc)='b' 

cbb=lc(chc); 

Cb(cdofl(chc),cdofl(chc))=Cb(cdofl(chc),cdofl(chc))  +  ebb; 
end 

changec=input('Are  there  any  other  changes  to  the  damping  matrix?  l=yes  2=no  '); 
end 
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changem=inputCWould  you  like  to  change  your  mass  matrix?  l=yes  2=no  '); 
while  changem==l 

chm=chm+l; 

flag_m=l; 

lm(chm)=input('input  value  of  added  mass  (lbs-secA2/in)  '); 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 

%  comparison  of  change  dof  with  the  c-set  matrix  and  formation  of 

%  new  c-set  matrix. 

% 

if  mdof (chm)~=cset 

cset=[cset  mdof(chm)]; 
end 

%  Classical  method  of  reformulatting  the  [M]  matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm))  +  lm(chm); 

changem=input('Are  there  any  other  changes  to  the  mass  matrix?  l=yes  2=no  '); 
end 

%  Handling  of  spring  added  to  ground 

% 

ground  =  find(cset==0); 

cset(ground)  =  Q; 


%  Formulation  of  iset  matrix  which  is  the  set  of  dofs  other  than 

%  those  in  the  cset  matrix  where  response  information  is  desired. 

% 

iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in?  '); 

%  Formation  of  zset  (c  U  b)  and  eset  (i  U  c  U  b)  matrices. 

% 

zset=[cset  bset];  eset=[iset  cset  bset]; 

%  Choosing  the  FRF  interested  in 
% 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

frf_index=find(eset=resp_dof); 

%  Formation  of  [del_Kc],  [del_Cc],  [del_Mc],  [del_Zc],  [del_Zb] 

%  &  [del_Z] 

% 

del_Kc=zerosGength(cset));  %  intializes  the  change  matrices 

del_Cc=zeros(length(cset));  %  to  zero  and  sets  their  sizes 

del_Mc=zeros(length(cset));  %  as  [cset  x  cset]  matrix 

ifflag_k=l 
fora=l:chk 
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if  (find(kdof l(a)=bset))~=[]  &  kdof2(a)='b' 

indexkl=[|;  %  provides  for  base  changes  to  not  be 

indexk2=Q;  %  included  in  del_Kc  matrix 

else 

indexk  1  =find(kdof  1  (a)=cset);    %  finds  where  in  cset  matrix 
indexk2=find(kdof2(a)=cset);    %  the  change  in  k  occurs 

end 

del_Kc(indexkl, indexk  l)=del_Kc(indexkl, indexk  1)  +  lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2)  -  lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl)  -  lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2)  +  lk(a); 
end 
end 

if  flag_c=l 
fora=l:chc 

if  (find(cdofl(chc)=bset))~=[]  &  cdof2(a)=='b* 
indexcl=[]; 
indexc2=[j; 
else 

indexcl=find(cdof  1  (a)=cset);    %  finds  where  in  cset  matrix 
indexc2=find(cdof2(a)==cset);    %  the  change  in  C  occurs 
end 

del_Cc(indexcl,indexcl)=del_Cc(indexcl,indexcl)  +  lc(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexcl,indexc2)  -  lc(a); 
del_Cc(indexc2,indexcl)=del_Cc(indexc2,indexcl)  -  lc(a); 
del_Cc(indexc2,indexc2)=del_Cc(indexc24ndexc2)  +  lc(a); 
end 
end 

if  flag_m==l 
fora=l:chm 
indexm=find(mdof(a)=cset);     %  finds  where  in  cset  matrix 

%  the  change  in  k  occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm)  +  lm(a); 
end 
end 

del_Kb=diag(kb);  %  diag  is  used  in  the  event  >1  base 

excitation 

del_Cb=diag(cb);  %  is  present 

% 

%  Time  Step: 

start_t  =  0.0; 
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%del_t   =0.0005;  %  plate  times 

%end_t    =.04; 

del_t   =0.0003;  %  beam  times 

end_t   =.03; 

%del_t   =  0.001;  %  spring  times 

%end_t   =.2; 

time  =  [start_t:del_t:end_t];    %  Time  points 
nstep  =  length(time);  %  No.  Time  points 

% 

%  Calculate  Impulse  Response  Functions  : 

% 


syn_t0  =  clock; 

[himp]  =  impmodal(wn,phi,zeta,Mmodal,time,zset); 

clear  wn  phi  Mmodal  zeta 

% 

%  Setup  and  Solve  Integral  Equation  for  x2*(t): 

% 

A  =  zeros(nstep); 

globalA  =  zeros(length(zset)*size(A,l)); 

count=0; 

col  =  l+count:nstep+count; 

for  acnt  =  1  :size(himp,  1 ); 

for  icnt_rows  =  2  :  nstep; 
[wts]  =  fTrapzWts(icnt_rows); 

A(icnt_rows,l:icnt_rows)  =  del_t  *  wts  .*  fliplr(himp(acnt,l:icnt_rows)); 
end 


row  =  col(l)4count:col  (length  (col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count  =  count  +nstep; 

if  row(length(row))  =  size(globalA) 

col  =  col  +  nstep; 

count  =  0; 
end 

end 
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clear  A 


%  Checks  to  see  if  there  are  any  redundant  dof  s  in  the  eset  due  to  changes  and 

%  bset  occurring  at  the  same  dofs.  If  there  is  redundancy,  it  is  located  in 

%  eset,  and  these  rows  and  columns  are  deleted  from  H_star.  This  is  done 

%  because  the  redundant  dof  column(s)  and  row(s)  in  H_star  creates  a  singular 

%  matrix  which  is  not  invertable  and  Kred  cannot  be  determined. 

% 

for  u=l  :length(bset) 

locate=find(bset(u)=zset); 

if  length(locate)  >  1 
pull=[pull  locate(2)]; 
redundant=[redundant  locate(l)]; 

end 
end 

pull_start=pull*nstep-(nstep- 1 ); 

pull_end=pull*nstep; 

delete_dof=[]; 

for  v=l:length(pull) 

delete_dof=[delete_dof  pull_start(v):pull_end(v)] ; 
end 

globalA(delete_dof,:)=[]; 
globalA(:,delete_dof)=Q ; 


% 


ff  =  ones(size(globalA,  1 ),  1); 

Yo  =  1 .0;  %  Amplitude  of  base  motion 

[Y,Ydot]  =  fBlastForcing(Yo,time',  'blst',  0); 

tol  =  le-4; 
dif=100; 

for  icnt=  1:300 

X  =  -globalA*ff; 

ificnt>=2 

[ii,jj]  =  max(abs(xlastl  -  X)); 

dif  =  max(abs(xlastl(jj)  -  X(jj))); 
end 

xlastl  =  X; 

if  dif  <  tol 

disp('Breaking');      break 
end 

[ff]  =  Synth_force(X,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t,bset,cset); 
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end 

syn_tl  =  clock; 
time_syn_time=etime(sy  n_t  1  ,syn_tO) 


% 


%  Classical  Method  to  Solve  for  Response: 

% 

clas_tO  =  clock; 

Xchk  =  zeros(length(zset),length(time)); 

[wn_c,phi_c,zeta_c,Mmodal_c,phi_norm_c]  =  modal(M_c,K_c,C_c); 

%zeta_c=0; 


for  icnt_modes  =  1  :  ndof; 

zeta_mode=zeta_c(icnt_modes); 

[mode_irf]  =  fModalIRF(wn_c(icnt_modes),  zeta_mode,  time); 
fmodal     =  phi_norm_c(:,icnt_modes)'  *  K_vec  *  Y; 
Xchk  =  Xchk  +  phi_norm_c(zset,icnt_modes)  *  fConvTrap(mode_iif,fmodaT,del_t); 
end 
clas_tl  =  clock; 
clas_syn_time=etime(clas_tl  ,clas_tO) 

%  Plotting  of  Transient  Response  using  Classical  &  Frequency  Synthesis  Methods 

% 

resp_index=find(zset==resp_dof); 

plot(time,Xchk(resp_index(l),:),,^~^time^C((l:nstep)+(^esp_index(l)-l)*nstep),,b,) 

grid 

title([Transient  Time  Response  at  dof  \int2str(resp_dof)]) 

ylabel('displacement  (in)') 

xlabel('time  (sec)') 

legendCClassicalVSynthesis') 
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tsynode45.m 


%  The  purpose  of  this  program  is  to  perform  the  time  domain  synthesis  on  a  structure  and 
%  compare  the  dynamic  responses  to  those  calculated  using  a  direct  integration  method. 

% 

clear 

disp('Select  the  file  which  contains  your  [M],  [K],  and  [C]  matrices' ) 
pause(2) 
[M_K_C,p]=uigetfile(,*.mat',,Load  [M],  [K],  [CD; 
load  (M_K_C) 
dofs=l:ndof; 
C=.00*C; 

[wn,phi,zeta,Mmodal,phi_norm]  =  modal(M,K,C); 

changek=2;   changec=2;     changem=2;     %  initializes  logic  variables 

%  to  'no'  status 

flag_k=0;    flag_c=0;     flag_m=0; 

chk=0;    chc=0;     chm=0;     %  initializes  counter  for  #  of  changes  to  matrices 


cset=0; 


%  initializes  cset  matrix  to  zero  to  avoid  null  index  error  in 
%  handling  the  spring  to  ground  if  no  changes  are  made  to  the 
%  structure. 


Hey_star=[J; 


HV_c=[]; 


K_c=K;    M_c=M;      C_c=C; 


K_vec  =  zeros(2*size(K_c,l),l); 
C_vec  =  zeros(2*size(C_c,l),l); 

Kb=zeros(size(K,  1 )); 
Cb=zeros(size(C,l)); 

b=num2strCb'); 


clearMKC 


%  (re)initializes  these  matrices  to  an  empty  matrix 
%  to  avoid  matrix  size  differences  each  time  the 
%  model  is  changed. 

%  initializes  change  matrices  which  will  be  used  in 
%  the  classical  method,  to  the  original  matrices. 


%  initializes  the  spring,  and  damper  to  base  matrices 
%  to  zero. 

%  allows  b  to  be  entered  as  a  numerical  value,  but 
%  declares  b  a  string  variable  to  used  for 
%  comparisons  in  logic  statements. 


%  Designation  where  base  excitation  is  located,  and  the  spring  constant  which 
%  connects  the  substructure  to  the  base  which  is  moving. 
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% 

bset=input('At  what  dof(s)  are  your  structure  excited?  '); 
for  b_spring=l:length(bset); 
kb(b_spring)=dnput(fprintf('input  the  value  of  the  spring  constant  which  connects  dof 
%g  to  the  base  ',bset(b_spring))); 
end 

cb=zeros(lJength(bset));  %  initializes  the  damper  constant  which  connects  the 

%  substructure  to  the  base  which  is  moving,  to  zero. 

%  Correcting  the  classical  [K]  &  [Kb]  matrices  to  include  the  spring  element 

%  connecting  the  substructure  to  the  base(s). 

% 

for  b_spring=l:length(bset); 

K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+kb(b_spring); 
Kb(Dset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring))  +  kb(b_spring); 
K_vec(size(K_c,  1  )+bset(b_spring))=K_vec(size(K_c,  1  )+bset(b_spring))  +  kb(b_spring); 
end 

changek=inputCWould  you  like  to  change  your  stiffness  matrix?  l=yes  2=no  '); 
while  changek=l 

chk=chk+ 1 ;  %  counter  for  determining  #  columns  in  mapping  matrix. 

flag_k=  1 ;  %  mapping  matrix  flag 

lk(chk)=input('input  value  of  added  spring  (lbs/in)  '); 

kdofl(chk)=input('input  constrained  dof  where  1st  end  of  added  spring  is  applied  '); 

kdof2(chk)=input('input  dof  where  2nd  end  of  added  spring  is  applied,  b  for  base,  or 

0  for  ground  '); 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 

%  are  not  included  in  the  bset.  A  change  from  dof  to  base  constitutes  an  addition  of 

%  excitation  points.  This  should  handled  when  asked  where  the  structure  is  excited. 

% 

if  kdofl(chk)~=bset  &  kdof2(chk)='b' 

while  kdofl(chk)~=bset  &  kdof2(chk)=='b* 
disp('Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
dispCIf  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change=l 

error('Go  back  and  change  your  basic  model.') 
else 

kdof  l(chk)=input('input  constrained  dof  where  1st  end  of  added  damper  is 

applied  '); 
kdof2(chk)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b 
for  base,  or  0  for  ground  '); 
end 
end 
end 

%  updates  the  base  spring  constant  and  does  not  count  these  changes 
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%  in  the  cset  matrix. 
% 

if  (find(kdofl(chk)==bset))~=D  &  kdof2(chk)=='b' 
base=find(kdof  1  (chk)==bset); 
kb(base)  =  kb(base)  +  lk(chk); 

K_vec(size(K_c,l)+kdofl(chk))=K_vec(size(K_c,l)+kdofl(chk))  +  lk(chk); 
else 
if  kdof  1  (chk)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  kdof  1  (chk)];  %  matrix  and  formation  of  new  c-set  matrix, 

end 

if  kdof2(chk)~=cset 
'     cset=[cset  kdof2(chk)]; 
end 
end 

%  Classical  method  of  reformulatting  the  [K]  matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk))  +  lk(chk); 

if  kdof2(chk)~=0  &  kdof2(chk)~=b' 

K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)Jcdof2(chk))  -  lk(chk); 

K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)Jcdofl(chk))  -  lk(chk); 

K_c(kdof2(chk)0cdof2(chk))=K_c(kdof2(chk)  Jcdof2(chk))  +  lk(chk); 
end 

%  Classical  method  of  reformulatting  [Kb]  matrix. 

% 

ifkdof2(chk)=*b' 

kbb=lk(chk); 

Kb(kdofl(chk),kdofl(chk))=Kb(kdofl(chk),kdofl(chk))  +  kbb; 
end 

changek=input('Are  there  any  other  changes  to  the  stiffness  matrix?  l=yes  2=no  '); 

end 

changec=input('Would  you  like  to  change  your  damping  matrix?  l=yes  2=no  '); 
while  changec=l 

chc=chc+l ;  %  counter  for  determining  #  columns  in  mapping  matrix. 

flag_c=l;        %  mapping  matrix  flag 

lc(chc)=input('input  value  of  added  damper  (lbs-sec/in)  '); 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is  applied  '); 
cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b  for  base,  or 

0  for  ground  '); 

%  Let's  the  user  know  that  he  is  not  allowed  to  make  dof(s)  to  the  base  changes  which 

%  are  not  included  in  the  bset.  A  change  from  dof  to  base  constitutes  an  addition  of 

%  excitation  points.  This  should  handled  when  asked  where  the  structure  is  excited. 

% 

if  cdofl(chc)~=bset  &  cdof2(chc)=='b' 

while  cdofl(chc)~=bset  &  cdof2(chc)=='b' 
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disp('Connecting  this  dof  to  the  base  is  changing  the  basic  model.') 
disp('If  this  is  what  you  desire,  go  back  and  include  this  dof  as  an  excitation 

point.') 
model_change=inputCWas  this  change  correct?  l=yes  2=no  '); 

if  model_change==l 

error('Go  back  and  change  your  basic  model.') 
else 

cdof  l(chc)=input('input  constrained  dof  where  1st  end  of  added  damper  is 
applied  '); 

cdof2(chc)=input('input  dof  where  2nd  end  of  added  damper  is  applied,  b 
for  base,  or  0  for  ground  '); 
end 
end 
end 

%  Updates  the  base  damper  constant  and  does  not  count  these  changes 

%  in  the  cset  matrix. 

% 

if  (find(cdofl(chc)=bset))~=[]  &  cdof2(chc)='b' 

base=find(cdof  1  (chc)=bset) ; 

cb(base)  =  cb(base)  +  lc(chc); 
C_vec(size(C_c,l)+cdofl(chc))=C_vec(size(C_c,l)+cdofl(chc))  +  lc(chc); 


else 


if  cdof  1  (chc)~=cset  %  comparison  of  change  dof  with  the  c-set 

cset=[cset  cdof  1  (chc)];       %  matrix  and  formation  of  new  c-set  matrix, 
end 

if  cdof2(chc)~=cset 
cset=[cset  cdof2(chc)]; 
end 


end 


%  Classical  method  of  reformulatting  the  [C]  matrix. 

% 

C_c(cdofl(chc),cdofl(chc))=C_c(cdofl(chc),cdofl(chc))  +  lc(chc); 

if  cdof2(chc)~=0  &  cdof2(chc)~='b' 

C_c(cdofl(chc),cdof2(chc))=C_c(cdofl(chc),cdof2(chc))  -  lc(chc); 

C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc))  -  lc(chc); 

C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc))  +  lc(chc); 
end 

%  Classical  method  of  reformulating  [Cb]  matrix. 

% 

ifcdof2(chc)=='b' 

cbb=lc(chc); 

Cb(cdofl(chc),cdofl(chc))=Cb(cdofl(chc),cdofl(chc))  +  ebb; 
end 

changec=input('Are  there  any  other  changes  to  the  damping  matrix?  l=yes  2=no  '); 
end 
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changem=inputCWould  you  like  to  change  your  mass  matrix?  l=yes  2=no  '); 
while  changem=l 

chm=chm+l; 

flag_m=l; 

lm(chm)=input('input  value  of  added  mass  (lbs-secA2/in)  '); 

mdof(chm)=input('input  constrained  dof  where  mass  is  added  '); 

%  comparison  of  change  dof  with  the  c-set  matrix  and  formation  of 

%  new  c-set  matrix. 

% 

if  mdof (chm)~=cset 

cset=[cset  mdof(chm)]; 
end 

%  Classical  method  of  reformulatting  the  [M]  matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm))  +  lm(chm); 

changem=input('Are  there  any  other  changes  to  the  mass  matrix?  l=yes  2=no  '); 
end 

%  Handling  of  spring  added  to  ground 

% 

ground  =  find(cset=0); 

cset(ground)  =  Q; 


%  Formulation  of  iset  matrix  which  is  the  set  of  dofs  other  than 

%  those  in  the  cset  matrix  where  response  information  is  desired. 

% 

iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in?  '); 

%  Formation  of  zset  (c  U  b)  and  eset  (i  U  c  U  b)  matrices. 

% 

zset=[cset  bset] ;  eset=[iset  cset  bset] ; 

%  Choosing  the  FRF  interested  in 
% 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

frf_index=find(eset==resp_dof); 

%  Formation  of  [del_Kc],  [del_Cc],  [del_Mc],  [del_Zc],  [del_Zb] 
%  &  [del_Z] 

% 

del_Kc=zeros(length(cset));  %  intializes  the  change  matrices 

del_Cc=zeros(length(cset));  %  to  zero  and  sets  their  sizes 

del_Mc=zerosQength(cset));  %  as  [cset  x  cset]  matrix 

ifflag_k==l 
fora=l:chk 

if  (find(kdofl(a)=bset))~=[]  &  kdof2(a)==*b' 
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indexkl=[];  %  provides  for  base  changes  to  not  be 

indexk2=[j;  %  included  in  del_Kc  matrix 

else 

indexkl=find(kdofl(a)=cset);    %  finds  where  in  cset  matrix 
indexk2=find(kdof2(a)=cset);    %  the  change  in  k  occurs 

end 

del_Kc(indexkl,indexkl)=del_Kc(indexkl,indexkl)  +  lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2)  -  lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl)  -  lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2)  +  lk(a); 
end 
end 

if  flag_c=l 
fora=l:chc 

if  (find(cdofl(chc)=bset))~=n  &  cdof2(a)=='b' 
indexcl=[]; 
indexc2=[]; 
else 

indexc  1  =find(cdof  1  (a)==cset);    %  finds  where  in  cset  matrix 
indexc2=find(cdof2(a)=cset);    %  the  change  in  C  occurs 
end 

del_Cc(indexc  1  ,indexc  1  )=del_Cc(indexc  1  ,indexc  1 )  +  lc(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexcl,indexc2)  -  lc(a); 
del_Cc(indexc2,indexcl)=del_Cc(indexc2,indexcl)  -  lc(a); 
del_Cc(indexc2,indexc2)=del_Cc(indexc24ndexc2)  +  lc(a); 
end 
end 

if  flag_m==l 
fora=l:chm 
indexm=find(mdof(a)==cset);     %  finds  where  in  cset  matrix 

%  the  change  in  k  occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm)  +  lm(a); 
end 
end 

del_Kb=diag(kb);  %  diag  is  used  in  the  event  >1  base  excitation 

del_Cb=diag(cb);  %  is  present 


% 


%  Time  Step: 

start_t  =  0.0; 

%del_t   =0.00025;  %  plate  times 

%end_t    =.02; 
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del_t   =0.0003;  %  beam  times 

end_t   =.03; 

%del_t   =  0.001;  %  spring  times 

%end_t   =  .2; 

time  =  [start_t:del_t:end_t];    %  Time  points 
nstep  =  length(time);  %  No.  Time  points 

% 

%  Calculate  Impulse  Response  Functions  : 

syn_t0  =  clock; 

[nimp]  =  impmodal(wn,phi,zeta,Mmodal,time,zset); 

clear  wn  phi  Mmodal  zeta 

% 

%  Setup  and  Solve  Integral  Equation  for  x2*(t): 

A  =  zeros(nstep); 

globalA  =  zeros(length(zset)*size(A,l)); 

count=0; 

col  =  l+count:nstep+count; 

for  acnt  =  1  :size(himp,l); 

for  icntjrows  =  2  :  nstep; 
[wts]  =  fTrapzWts(icnt_rows); 

A(icnt_rows,l:icnt_rows)  =  del_t  *  wts  .*  fliplr(himp(acnt,l:icnt_rows)); 
end 


row  =  col(l)+count:colGength(col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count  =  count  +nstep; 

if  row(length(row))  ==  size(globalA) 

col  =  col  +  nstep; 

count  =  0; 
end 


end 

clear  A 

%disp(sprintf('Norm(globalA)  =  %5.3f  ,norm(globalA))) 
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% 

ff  =  ones(length(zset)*nstep,l); 

Yo  =  1 .0;  %  Amplitude  of  base  motion 

[Y,Ydot]  =  fBlastForcing(Yo,time\  'blst',  0); 

tol  =  le-4; 
dif=100; 

for  icnt=  1:300 

X  =  -globalA*ff; 

if  icnt>=2 

[ii,jj]  =  max(abs(xlastl  -  X)); 

dif  =  max(abs(xlastl(jj)  -  XQj))); 
end 

xlastl  =  X; 

if  dif  <  tol 

disp('Breaking');      break 
end 


[ff]  =  Synth_force(X,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t,bset,cset); 


end 


syn_tl  =  clock; 
time_syn_time=etime(syn_t  1  ,syn_t0) 


%  Classical  Method  to  Solve  for  Response: 

clas_t0  =  clock; 

save  structure  K_c  C_c  M_c  K_vec  C_vec  Yo 

Xchk0=zeros(2*ndof ,  1 ); 

[tchk,Xchk]=ode45(,scnd_to_frst,,start_t,end_t,XchkO); 

clas_tl  =  clock; 

clas_sy  n_time=etime(clas_t  1  ,clas_t0) 


%  Plotting  of  Transient  Response  using  Classical  &  Frequency  Synthesis  Methods 

% 

resp_index=find(zset=resp_dof); 
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plot(tchk,Xchk(:,resp_dof(l)),,r--,,time,X((l:nstep)+(resp_index(l)-l)*nstep),'b') 

grid 

title([Transient  Time  Response  at  dof  ',int2str(resp_dof)]) 

ylabel('displacement  (in)') 

xlabel('time  (sec)') 

legend('Classical','Synthesis') 
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tsynmain.m 


%  This  program  is  designed  to  calculate  the  new  transient  responses  of  a  system  after 

%  changes  in  the  mass,  stiffness,  and  damping  matrices  have  been  made. 

%  The  new  responses  will  be  calculated  using  the  synthesis  method  in  the  time  domain. 

% 

clear 

load  change 

%  Designation  of  the  frequency  range  and  step  size 

% 

%  initial      step  size      end 

t_input  =  [  0.0  0.0004  .07  ]; 

time  =  t_input(l):t_input(2):t_input(3); 

nstep  =  length(time);  %  No.  Time  points 

[wn,phi,zeta,Mmodal,phi_norm]  =  modal(M,K,C); 
clear  M  K  C 

syn_t0  =  clock; 

[himp]  =  impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear  wn  phi  zeta  Mmodal 

[globalA]  =  buildA(himp,bset,zset,nstep,t_input(2)); 

Yo  =  1.0;  %  Amplitude  of  forcing  function 

[Y,Ydot]  =  fBlastForcing(Yo,time*,  'blst',  0); 

[X_star,icnt]  =  timesynth(globalA,zset,Y,Ydot,del_Kc,del_Cc,del_McJcb,cb,... 

nstep,t_input(2),bset,cset); 

syn_tl  =  clock; 
time_syn_time=etime(syn_t  1  ,syn_t0) 

%  Plotting  of  Transient  Response  using  Classical  &  Frequency  Synthesis  Methods 

% 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

[X_star_desired]  =  time_sift(zset,resp_dof,X_star,nstep); 

plot(time,X_star_desired,'b') 

grid 

title( ['Synthesized  Transient  Time  Response  at  dof  ',int2str(resp_dof)]) 

ylabel('displacement  (in)') 
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xlabel('time  (sec)') 
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impmodal.m 


function  [himp]  =  impmodal(wn,phi,zeta,Mmodal,t,imp_dof) 

%  This  function  is  designed  to  calculate  the  IRF  using  mode  shapes  and 
%  summing  the  IRF  over  a  number  of  modes. 

% 

ncol=(length(imp_dof)*  (tength(imp_dof)+ 1  ))/2; 

%  Calculation  of  impulse  response  function 

% 

himp=zeros(ncol,length(t)); 


nmodes=size(phi);       b=0; 


forb=l:nmodes 


%  Initialization  to  zero,  and  sizing 
%  of  the  freq  x  lower  triangle  matrix. 

%  Defining  the  number  of  modes  that  exist  and 
%  initializing  b  which  keeps  track  of  the 
%  of  modes 


%  Elimination  of  unnessecary  elements  and  rearrangement  of  original 

%  mode  shapes  {phi}  to  form  {phi_red}  and  [num_red]  which  is  now  an 

%  imp_dof  x  imp_dof  matrix. 

% 

num_red=phi(imp_dof,b)*phi(imp_dof,b)'; 


ifwn(b)>le-3 


%  Then  elastic  mode 


wd  =  wn(b)  *  sqrt(l  -  zeta(b)A2); 

mode_irf  =  exp(-zeta(b)*wn(b)*t)  .*  sin(wd  *  t)  /  (Mmodal(b)*wd); 


%  Rigid  body  mode 


elseif  wn(b)  <=  le-3 
mode_irf  =  t/Mmodal(b); 
end 


%  Changes  the  [num_red]  matrix  from  an  imp_dof  x  imp_dof  matrix  to 
%  a  1  x  lower  triangle  vector. 


% 
H_lowtri=tril(num_red); 

symtry  =  find(H_lowtri==0); 

H_lowtri(symtry)  =  []; 


%  pulls  out  the  lower  triangle 

%  and  places  zeros  everywhere  else 

%  finds  zero  values  in  the 
%  lower  triangle  matrix 

%  deletes  all  values  of  zero  and 


175 


%  turns  the  remaining  elements 
%  into  a  1  x  lengthOower  triangle) 
%  vector 


for  cnt=l  :length(H_lowtri) 

himp_mode(cnt,:)=H_lowtri(cnt)*mode_irf; 
end 

himp  =  himp  +  himp_mode; 
end 
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buildA.m 


function  [globalA]  =  buildA(himp,bset,zset,nstep,del_t) 

%  The  purpose  of  this  program  is  to  build  the  trapezoidal  quadrature  matrices  to  be  used 
%  for  the  solution  of  the  VIDE 

% 

A  =  zeros(nstep); 

globalA  =  zeros(length(zset)*size(A,l)); 

count=0; 

col  =  l+count:nstep+count; 

for  acnt  =  l:size(Mmp,l); 

for  icnt_rows  =  2  :  nstep; 
[wts]  =  fTrapzWts(icnt_rows); 

A(icnt_rows,l:icnt_rows)  =  del_t  *  wts  .*  fliplr(himp(acnt,I:icnt_rows)); 
end 

row  =  col(l)+count:colGength(col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count  =  count  +nstep; 

if  row(length(row))  =  size(globalA) 
col  =  col  +  nstep; 
count  =  0; 
end 
end 

%  Checks  to  see  if  there  are  any  redundant  dofs  in  the  eset  due  to  changes  and 

%  bset  occurring  at  the  same  dofs.  If  there  is  redundancy,  it  is  located  in 

%  eset,  and  these  rows  and  columns  are  deleted  from  H_star. 

% 

for  u=l  :length(bset) 

locate=find(bset(u)=zset); 

if  length(locate)  >  1 
pull=[pull  locate(2)]; 
redundant=[redundant  locate(l)]; 

end 
end 

pull_start=pull*nstep-(nstep- 1 ); 

pull_end=pull*nstep; 

delete_dof=Q; 

for  v=l:length(pull) 

delete_dof=[delete_dof  pull_start(v):pull_end(v)] ; 
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end 

globalA(delete_dof,:)=[]; 
globalA(:,delete_dof)=Q; 
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fBlastForcing.m 


function  [f_of_t,fdot]  =  fBlastForcing(Fo,time,  type,  plotit); 

% 

%    Usage:  [f_of_t,fdot]  =  fBlastForcing(Fo,time,  type,  plotit); 

% 

%  Choices:  sine   blst    step 

% 

%      type  =  'step'    STRING  Variable 

% 

% 

%     If  use  'sine1,  fdot  also  returned. 

% 

%    This  function  returns  a  forcing  function  which  is 

%     a  "blast"  function. 

% 

%         F(t)  =  Fo  *  (  exp(-at)  -  exp(-bt) ) 

% 

%    where  a  and  b  are  constants  which  shape  the  blast, 

%    and  Fo  is  the  amplitude  of  the  blast. 

% 

%    The  variable  "plotit"  is  a  switch  which  if  set  =  1  will 

%    cause  the  f(t)  to  be  plotted,  if  set  to  anything  else 

%     will  not  plot 

% 

% 

%  Choices:  sine   blst    step 

%    type  =  'step'; 

if  type  =  'blst'; 

%       dispC  Blast  forcing  used...') 
a  =  100.0; 
b  =  300.0; 

f_of_t  =  Fo  *  ( exp(-a*time)  -  exp(-b*time) ); 
fdot  =  Fo  *  ( -a*exp(-a*time)  +  b*exp(-b*time) ); 

elseif  type  ==  'step'; 

%       dispC  Step  forcing  used...') 
f_of_t  =  Fo  *  ones(size(time)); 
fdot  =  []; 

elseif  type  ==  'sine'; 
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%       dispC  Sine  forcing  used...') 
W  =  5;  %  Hertz 

f_of_t  =  Fo  *  sin(2*pi*W*time); 
fdot   =  Fo  *  (2*pi*W)*cos(2*pi*W*time); 

end; 

if  plotit  =  1; 

figure(gcf+l) 

if  type  =  'sine'; 

plot(time,f_of_t,time,fdot);grid 
else 
plot(time,f_of_t);grid 
end 
pause 
%    elf 
end 

%  End  function. 
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timesynth.m 


function  [X_star,icnt]  =  timesynth(globalA,zset,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,. 

nstep  ,dd_t,bset,cset) 

%  This  function  is  designed  to  calculate  the  new  dynamic  response  using  the  synthesis 

%  method. 

% 

ff  =  ones(size(globalA,l),l); 

tol  =  le-2; 

dif=100; 

for  icnt=  1:300 

X_star  =  -globalA*ff; 

if  icnt  >=  2 

[iijj]  =  max(abs(xlastl  -  X_star)); 

dif  =  max(abs(xlastl(jj)  -  X_star(jjj))); 
end 

ificnt>=300 
ii 

jj 
dif 
end 

xlastl  =  X_star; 

if  dif  <  tol 
%  ciisp(,B^eaking,); 

break 

end 
%         icnt 

[ff]  = 
Synth_force(X_star,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t,bset,cset); 

end 
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time  sift.m 


function  [X_star_desired]  =  time_sift(zset,resp_dof,X_star,nstep) 

%  The  purpose  of  this  program  is  to  sort  through  the  synthesized  transient 
%  responses  and  return  the  one  in  which  the  user  is  interested  in. 

% 

resp_index=find(zset=resp_dof); 

X_star_desired  =  X_star((l:nstep)+(resp_index(l)-l)*nstep); 
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APPENDIX  D.  OPTIMIZATION  COMPUTER  CODES 

The  following  is  a  list  and  a  brief  description  of  the  main  MATLAB 
computer  codes  that  were  written  in  order  to  perform  the  optimization  of  a 
shock  and  vibration  isolation  system. 

•  isomain.m  -  the  main  program  which  loads  the  baseline  structure  to  be 
optimized,  and  allows  for  the  designation  of  the  design  variables.  Also,  if 
necessary,  calls  various  modularized  functions  such  as  modal.m, 
frfmodal.m,  impmodal.m,  fBlastForcing.m,  buildA.m.  It  then  allows  for 
the  optimization  inputs,  the  use  of  the  MATLAB  Optimization  function 
constr.m,  and  the  post  processing. 

•  modal.m  -  solves  for  the  structure's  natural  frequencies,  mode  shapes  and 
other  modal  matrices  and  vectors. 

•  frfmodal.m  -  calculates  the  structure's  FRFs. 

•  impmodal.m  -  calculates  the  structure's  IRFs. 

•  buildA.m  -  builds  the  quadrature  matrices. 

•  fBlastForcing.m  -  inputs  the  base  excitation  as  a  function  of  time. 

•  isosub.m  -  the  subroutine  program  which  is  called  by  the  MATLAB 
Optimization  function  constr.m.  This  is  where  the  optimization  iterations 
occur,  the  design  variables  change,  and  the  system's  responses  are 
calculated.  This  program  calls  the  functions  delta.m,  frfsynth.m, 
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statdelta.m,  statsynth.m  and  timesynth.m.  It  then  calculates  the 
normalized  constraint  values. 

•  delta.m  -  forms  the  modification  matrices  which  will  be  used  to  form 
impedance  matrices  or  determine  coupling  forces. 

•  frfsynth.m  -  performs  the  synthesis  on  the  baseline  structure  and  returns 
the  synthesized  FRFs. 

•  statdelta.m  -  forms  the  modification  matrices  which  will  be  used  to  form 
impedance  matrices. 

•  statsynth.m  -  performs  the  synthesis  on  the  baseline  structure  and  returns 
the  synthesized  static  displacements. 

•  timesynth.m  -  performs  the  synthesis  on  the  baseline  structure  and 

returns  the  synthesized  transient  responses. 

The  full  codes  are  contained  on  subsequent  pages.  Any  codes  that  were 
previously  presented  will  not  be  repeated. 
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isomain.m 


clear 

clear  global 
diary  isomain.dia 
tic 


global  kdofl  kdof2  cdofl  cdof2  mdof  iset  cset  bset  zset  eset  kbO  cbO  hee  omega ... 
opt_omega  excite  resp_dof  inp_dof  resp_index  stat_omega  stat_hee  iter ... 
del_vars  del_f  H_desired_0  del_H_desired  globalA  Y  Ydot  t_input  nstep 

% 

%         Loading  of  Structure  to  be  Optimized 

% 

disp('Are  there  already  reduced  FRF  &  IRF  matrices  in  the  correct  format  available,') 

disp('or  does  the  presynthesized  FRF  &  IRF  matrices  need  to  be  generated  using  the  dof ') 

disp('locations  where  you  desire  to  change  the  structure?  ') 

FRFJRF=input('l  reduced  FRF/IRF  already  exists         2=generate  reduced  FRF/ERF  '); 

%  Protects  against  user  making  an  error  in  choosing  how  FRF/IRF  is  obtained. 
% 

ifFRF_IRF~=[12] 
whileFRF_IRF~=[12] 

disp('Error  in  choosing  how  FRF/IRF  is  obtained.  Choose  1  or  2.') 
FRF_IRF=input(*l=reduced  FRF/IRF  already  exists  2=generate  reduced  FRF/IRF  '); 
end 
end 

ifFRF_IRF=l 

disp('Select  the  file  which  contains  your  presynthesized  FRF/IRFs,  a  corresponding 
frequency  ') 

disp('and  time  vectors,  and  the  vector  of  all  the  dofs  that  will  be  in  your  eset  vector.') 

pause(2) 

[hee_dofs_omega,p]=uigetfile(,*.maf,'Load  FRF/IRF') 

load  (hee_dofs_omega) 

disp('Also  select  the  file  which  contains  your  presynthesized  FRFs  with  zero  damping  ') 

disp('for  the  purpose  of  calculating  static  displacement.') 

pause(2) 

[stat_hee,p]=uigetfile('*.mat','Load  Static  FRF*) 

load  (stat_hee) 

else 
disp('Select  the  file  which  contains  your  [M],  [K],  and  [C]  matrices' ) 
pause(2) 
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[M_K_C,p]=uigetfile(,*.mat\,Load  [M],  [K],  [C]'); 
load  (M_K_C) 
dofs=l:ndof; 
kcO=K(109,109); 
end 

% 

%         Designation  of  Optimization  and  Synthesis  DOFs 

% 

kdofl=[];  kdof2=[];         %  initalizes  the  matrices  which  keep 

cdof  1 =D ;  cdof2=[] ;         %  track  of  the  dofs  where  changes  occur 

mdof=D; 

b=num2str('b');  %  allows  b  to  be  entered  as  a  numerical  value  for  the 

dof 

%  choices;  but  declares  b  a  string  variable  to 
used  for 

%  comparisons  in  logic  statements. 

kdofs=input('Input  dof  pairs  for  stiffness  changes  for  optimization  '); 
ifkdofs=[] 

kdofl=[];      kdof2=n; 
else 

kdofl=kdofs(:,l);  kdof2=kdofs(:,2); 

end 

cdofs=input(lnput  dof  pairs  for  damping  changes  for  optimization  '); 
if  cdofs=[] 

cdofl=[];      cdof2=D; 
else 

cdof  1  =cdofs(:,  1 );  cdof2=cdofs(:,2); 

end 

mdof=input('Input  dofs  for  mass  changes  for  optimization  '); 

% 

%         Formation  of  Partitioning  and  Arrangement  Vectors 

% 

% 

excite=inputCHow  is  the  structure  excited?  l=system    2=base  '); 

%  Formation  of  {bset} 

% 

cset=0;  bset=[];  %  initiahzes  cset  matrix  to  zero  to  avoid 

%  null  index  error  if  no  changes  are  made  to 
%  the  structure  and  bset  to  empty  matrix  to 
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%  prevent  error  in  the  event  this  is  not  a 
%  base  excitation  and  bset  does  not  exist 

kbO=0;  cbO=0;  %  initializes  kb  and  cb  to  zero  to  prevent  error  in 

%  the  event  this  is  not  a  base  excitation  and  kb  &  cb  does 
%  not  exist. 

if  excite=2 
%  Designation  where  base  excitation  is  located. 
% 
bset=input('At  what  dof(s)  are  your  structure  excited?  '); 

%  Designation  of  the  spring  and  damper  constants  which  connects  the  substructure  to  the 

%  base  which  is  moving.and  protection  against  user  making  an  error  in  choosing  the 

%  value  of  the  constant(s). 

% 

for  b_spring=l:length(bset) 

kbO(b_spring)=input(fprintf(,input  the  value  of  the  spring  constant  which  connects 

dof  %g  to  the  base  ',bset(b_spring)»; 
kbO_zero=find(kbO=0); 
if  length(kbO)~=b_spring  I  kbO_zero~=D 

while  length(kbO)~=b_spring  I  kbO_zero~=[] 

fprintf('You  made  an  error  in  entering  the  value  for  dof  %g  base  excitation 

spring  constant.\n',bset(b_spring)) 
kbO(b_spring)=0; 
kbO(b_spring)=input(fprintf('Re-input  the  value  of  the  spring  constant  which 

connects  dof  %g  to  the  base  \bset(b_spring))); 
kbO_zero=find(kbO==0); 
end 
end 
end 

cbO=zeros(  1  ,length(bset));  %  initializes  the  damper  constant  which  connects  the 

%  substructure  to  the  base  which  is  moving,  to  zero, 
end 

%  Formation  of  {cset}  vector 
% 

ifkdofs~=D 
for  chk=l:length(kdofl) 

%  Comparison  of  change  dof  with  the  c-set  matrix  and  formation  of  new  c-set  matrix. 

%  Does  not  count  changes  to  the  base  spring  constant  in  the  cset  matrix. 

% 

if  excite==2  &  (find(kdofl(chk)=bset))~=[]  &  kdof2(chk)=*b' 

cset=cset; 
else 

if  kdof  l(chk)~=cset 
cset=[cset  kdofl(chk)]; 

end 

if  kdof2(chk)~=cset 
cset=[cset  kdof2(chk)]; 
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end 
end 
end 
end 

ifcdofs~=[] 
for  chc=l  :length(cdof  1) 

%  Comparison  of  change  dof  with  the  c-set  matrix  and  formation  of  new  c-set  matrix. 
%  Does  not  count  changes  to  the  base  damper  constant  in  the  cset  matrix. 
% 

if  excite=2  &  (find(cdofl(chc)=bset))~=[]  &  cdof2(chc)='b' 
cset=cset; 


else 


if  cdof  1  (chc)~=cset 
cset=[cset  cdofl(chc)]; 

end 

if  cdof2(chc)~=cset 
cset=[cset  cdof2(chc)]; 

end 


end 
end 
end 


ifmdof~=[] 
for  chm=l:length(mdof) 

%  Comparison  of  change  dof  with  the  c-set  matrix  and  formation  of  new  c-set  matrix. 

% 

if  mdof(chm)~=cset 

cset=[cset  mdof(chm)];  %  and  formation  of  new  c-set  matrix, 

end 
end 
end 

ground  =  find(cset=0);  %  Handling  of  spring  added  to  ground  to 

eliminate 

cset(ground)  =  [];  %  zero  from  cset 

%  Formation  of  {iset}  vector 

% 

ifFRF_IRF=l 

disp('Since  you  inputed  your  own  FRF/IRF  matrix,  your  iset  is  chosen  as  all  dofs 
remaining ') 

disp('in  your  eset/dof  vector  after  extracting  the  cset  vector.*) 

iset=dofs; 

for  a=l  ilength(cset) 
extract(a)=find(cset(a)==dofs); 

end 

iset(extract)=[]; 
else 
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%   iset=inputCWhat  dofs  other  than  those  where  change  occured,  are  you  interested  in? 

*); 

iset=[]; 
end 

%  Formation  of  zset  (c  U  b)  and  eset  (i  U  c  U  b)  matrices. 

% 

zset=[cset  bset];  eset=[iset  cset  bset]; 

% 

%         Formulation  of  original  [H],  [Himp],  for  Freq  &  Time  Synthesis  and  Static 

%         Displacement 

% 


if  FRF_IRF=2 

%  Designation  of  the  frequency  and  time  range  and  step  size 

% 

%  initial      step  size      end 

freq  =  [      .01  10        le3+.01  ]; 

t_input  =  [  0.0  0.001  .05  ]; 

omega=rreq(l):freq(2):freq(3); 

time  =  t_input(l):t_input(2):t_input(3); 

nstep  =  length(time);  %  No.  Time  points 

[wn,phi,zeta,Mmodal,phi_norm,Cmodal]  =  modal(M,K,C); 
[nee]  =  frfmodal(wn,phi,zeta,Mmodal,omega,eset); 
[himp]  =  impmodal(wn,phi,zeta,Mmodal,time,zset); 

%  Information  for  static  displacement  synthesis 

% 
stat_zeta  =  0*zeta; 

stat_omega  =  .  1 : .  1 :  .4; 
[stat_hee]  =  frfmodal(wn,phi,stat_zeta,Mmodal,stat_omega,eset); 

end 

clear  M  K  C 

% 

%         Time  Synthesis 

%         Builds  the  matrix  A  which  uses  trapezoidal  weights  to  integrate  the  IRFs 

% 

[globalA]  =  buildA(himp,bset,zset,nstep,t_input(2)); 

clear  himp 
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%         Forms  the  base  excitation  vectors  of  displacement  and  velocity 

% 

Yo  =  1 .0;  %  Amplitude  of  forcing  function 

[Y,Ydot]  =  fBlastForcing(Yo,time',  'blst',  0); 

% 

%         The  Optimization 

% 


%         Designates  the  frequency  over  which  the  optimization  will  be  conducted 

% 

ifFRF_IRF=2 

disp('Since  you  generated  the  FRF/IRFs  through  this  optimization  program,  your  FRF') 

disp('frequencies  and  IRF  times  are  automatically  chosen  as  the  freq.  and  time  range  ') 

disp('for  your  optimization.') 

opt_omega=omega; 
else 

opt_freql=inputCWhat  is  the  initial  frequency  over  which  to  evaluate  this  problem'); 

opt_freq2=input('What  is  the  final  frequency  over  which  to  evaluate  this  problem1); 

disp(The  frequency  increment  is  automatically  chosen  to  match  the  FRF  freq. 
increment') 

opt_omega=opt_freq  1  :omega(2)-omega(  1  ):opt_freq2; 

end 

%         Selecting  the  FRF  that  is  to  be  minimized. 

% 

resp_dof=inputCWhat  response  dof  are  you  interested  in?  '); 

if  excite=2 

inp_dof=D; 
else 

inp_dof=inputCWhat  input  dof  are  you  interested  in?  '); 
end 
resp_index=find(zset=resp_dof); 

%         Setting  the  number  of  of  optimization  variables,  their  initial  values,  and 

%         their  upper  and  lower  bounds. 

% 

numoptvar=inputCHow  many  optimization  variables  do  you  have?  '); 

del0=[l  1  1  1]; 

%  del_kp  del_kc     del_mp 

% 

%  del_kb  del_kc  del_cb     del_cc 

dellb=[  -4  -4  0  0         ]; 

delub=[  5  5  5  5         ]; 
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%         Initializing  the  matirices  which  will  keep  track  of  how  the  design  variables 

%         and  objective  function  changes  for  every  optimization  iteration  for  the 

%         purpose  of  plotting. 

% 

del_vars  =  Q;   del_f  =  [];  del_H_desired  =  Q; 

iter=0;  %  Sets  the  number  of  iterations  counter  at  zero 

%         Calling  of  constr.m  routine  to  perform  optimization 
%         on  subroutine  iso_sub.m 

% 

options(l)=l; 

options(14)=200; 

toe 

tic 
[del]=constr('isosub,,delO,options,dellb,delub); 

%  OUTPUT 

load  opt_data 

final_del_kb  =  del_vars(length(del_vars),l); 
final_del_kc  =  del_vars(length(del_vars),2); 
final_del_cb  =  del_vars(length(del_vars),3); 
final_del_cc  =  del_vars(length(del_vars),4); 

disp(The  recommended  change  in  base  isolator  stiffness  (lb/in)  is  — >  '),final_del_kb 

disp(The  recommended  change  in  computer  isolator  stiffness  (lb/in)  is  — >  '),final_del_kc 

disp(The  recommended  change  in  base  isolator  damping  (lb/in/s)  is  — >  '),final_del_cb 

disp(The  recommended  change  in  computer  isolator  damping  (lb/in/s)  is  — > 

'),final_del_cc 

disp('Constraints:') 

disp(g) 

iterations=l  liter, 

%  Plot  of  FRF  to  observe  how  it  changes  during  the  optimization 

% 

figure(l) 

semilogy(omega,abs(del_H_desired(l,:)),,--r,,omega,abs(del_H_desired(2,:)),,m') 

title('Optimized  Frequency  Response  Function  [H*]') 

ylabel(TRF  Amplitude  (unitless)') 

xlabel(Trequency  (rad/s)') 

legendCBefore'/After') 

figure(2) 
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subplot(2,l,l) 

plot(iterations,del_vars(:,l),,-r',iterations,del_vars(:,2),,m') 

title('Change  in  Variables  vs.  Iterations') 

ylabelCIsolator  Stiffness  (lb/in)') 

legendCdeLkbVdeLkc') 
subplot(2,l,2) 

plot(iterations4el_vars(:3)/"r',iterations,del_vars(:,4),'m') 

title('Change  in  Variables  vs.  Iterations') 

xlabel('#  of  iterations') 

ylabel('Isolator  Damping  (lb-s/in)') 

legendCdeLcb'.'deLcc') 

figure(3) 
subplot(2,l,l) 

plot(iterations,del_vars(:,  1  )+kbO(  1  ),'--r',iterations,del_vars(:  ,2)+kc0,'m') 

title('Change  in  Isolation  Constants  vs.  Iterations') 

ylabel('Isolator  Stiffness  (lb/in)') 

legend('k_base','k_computer') 
subplot(2,l,2) 

plot(iterations,del_vars(:,3),,~r',iterations,del_vars(:,4),,m') 

title('Change  in  Isolation  Constants  vs.  Iterations') 

xlabel('#  of  iterations') 

ylabelCIsolator  Damping  (lb-s/in)') 

legend('c_base','c_computer') 

figure(4) 

semilogy(iterations,del_f,'r') 

title('Change  in  Objective  Function  vs.  Iterations') 

ylabel('Magnitude') 

xlabel('#  of  iterations') 

figure(5) 

plot(time,X_star((  1  :nstep)+(resp_index(l )- 1  )*nstep),'b') 

grid 

title(['Synthesized  Transient  Time  Response  at  dof  \int2str(resp_dof)]) 

ylabel('displacement  (in)') 

xlabel('time  (sec)') 

figure(6) 

plot(time,Xaccel_star/386.4,'m') 

grid 

title(['Synthesized  Transient  Time  Response  at  dof  \int2str(resp_dof)]) 

ylabel('acceleration  (g"s)') 

xlabel('time  (sec)') 

toe 

diary  off 

%  end  Iso  main.m 
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isosub.m 


function  [f,g]=iso_sub(del) 

global  kdofl  kdof2  cdofl  cdof2  mdof  iset  cset  bset  zset  eset  kbO  cbO  hee  omega... 

opt_omega  excite  resp_dof  inp_dof  resp_index  stat_omega  stat_hee  iter  del_vars... 
del_f  H_desired_0  del_H_desired  globalA  Y  Ydot  t_input  nstep 

iter=iter  +  1;  %  counts  number  of  iterations 

%  Initialization  and  assigning  of  the  optimization  variables  to  changes  which  are  made  to 

%  the  structure.  Spring  and  damper  constants  are  scaled  to  even  out  the  domain. 

% 

lk=0;    lc=0;     lm=0; 

del_kb=le3*del(l);  del_kc=le3*del(2); 

del_cb=l*del(3);  del_cc=l*del(4); 

lk(l  :4)=del_kb*ones(l  ,4);  lk(5)=del_kc; 

lk(6:9)=-le3*ones(l,4); 

lc(  1 :4)=del_cb*ones(  1 ,4) ;  lc(5)=del_cc; 

a **** ******* ************* *** ************* ** ***************** *********** 
% 

%         Formulation  of  FRF  and  Calculation  of  mean  square  acceleration 

[del_Kc,del_Cc,del_McJcb,cb]  = 

delta(cset,bset,kdof  1  ,kdof2,lk,cdof  1  ,cdof2,lc,mdof ,lm,kbO,cbO); 

[h_star]  =  frfsynth(hee,kb,cb,omega,excite,eset,iset,bset,del_Kc,deLCc,del_Mc); 

[H_desired]  =  frf_sift(eset,resp_dof,inp_dof,excite,h_star); 

cut_freq=[find(omega<opt_omega(  1 ))  find(omega>opt_omega(length(opt_omega)))] ; 
H_desired(cut_freq)=[] ; 
H_objective=(abs(H_desired)).A2; 
[f]=simpl_3(H_objective,opt_omega(2)-opt_omega(l)) 

en l  ********************************************************************** 
% 

%         Plot  Generation 

%  Saving  a  matrix  of  the  optimization  variables,  objective  function,  and  FRF  in  order  to 
%  plot  how  they  changed  during  the  optimizatiion. 

% 
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if  iter  =1 

del_vars  =  [del_kb  del_kc  del_cb  del_cc]; 

H_desired_0  =  H_desired; 

del_f=[f]; 
else 

del_vars  =  [del_vars;del_kb  del_kc  del_cb  del_cc]; 

del_H_desired  =  [H_desired_0;  H_desired]; 

del_f=[del_f;f]; 
end 

clear  hee  H_desired_0  H_desired 

% 

%         Calculation  of  plate  and  computer  static  defection 


[stat_del_Kc,stat_del_Mc,stat_kb]  =  statdelta(cset,bset,kdof  1  Jcdof2,lk,mdof,lmjcb0); 

[z_stat,Fred,Kred,Mred,HstarO_dp]  =  statsynth(stat_hee,stat_kb,stat_omega  ,... 

eset,iset,bset,stat_del_Kc,stat_del_Mc); 

% 

%  Calculation  of  dynamic  response  due  to  shock 

[X_star,icnt]  =  timesynm(globalA,zset,Y,Ydot,dd_Kc,dd_Cc,del_Mc,kb,cb,... 

nstep,t_input(2),bset,cset); 

fprintf('icnt  =  %g',icnt) 

%dispCFinish  time  synthesis') 

% 

%  Calculation  of  normalized  constraint  values 

maxzp_stat=-.08;  maxzc_stat=-.03; 

maxkp_stretch=l;  maxkc_stretch=l; 

maxzc_accel=30*  3  86.4; 

zp_stat=z_stat([l  3:6]);  %  pulls  out  the  static  displacements  for  the  plate 

zp_stat=-(max(abs(zp_stat)));  %  and  determines  where  the  maximum  displacement 

%  is  min  command  used  because  static  displacement 

%  is  negative 

zc_stat=-(abs(z_stat(2))  -  abs(z_stat(l)));  %  calculation  of  the  relative  static 

%  displacement  between  the  plate  and 
%  computer 
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%  calculation  of  the  maximum  relative  dynamic  displacements  between  the  plate  and  the 

%  base  and  the  computer  and  the  plate 

% 

kp_stretch  =  X_star(2*nstep+l:6*nstep)  -  [Y;Y;Y;Y]; 

kp_stretch  =  max(abs(kp_stretch)); 

kc_stretch  =  X_star(nstep+l:2*nstep)  -  X_star(l:nstep); 

kc_stretch  =  max(abs(kc_stretch)); 

%  calculation  of  the  maximum  absolute  acceleration  of  the  computer 

% 

Xaccel_star  =  fddot(X_star((  1  :nstep)+(resp_index(  1 )- 1  )*nstep),t_input(2)); 

zc_accel  =  max(abs(Xaccel_star)); 

g(l)=zp_stat/maxzp_stat  - 1; 
g(2)=zc_stat/maxzc_stat  - 1; 
g(3)=kp_stretch/maxkp_stretch  - 1; 
g(4)=kc_stretch/maxkc_stretch  - 1; 
g(5)=zc_accel/maxzc_accel- 1 ; 

save  opt_data  del_vars  del_f  iter  del_H_desired  g  z_stat  kp_stretch  kc_stretch  zc_accel . 
X  starXaccel  star 
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